Abstract:
A method solves a stochastic quadratic program (StQP) for a convex set with a set of general linear equalities and inequalities by an alternating direction method of multipliers (ADMM). The method determines an optimal solution, or certifies that no solution exists. The method optimizes a step size β for the ADMM. The method is accelerated using a conjugate gradient (CG) method. The StMPC problem is decomposed into two blocks. The first block corresponds to an equality constrained QP, and the second block corresponds to a projection onto the StMPC inequalities and anticipativity constraints. The StMPC problem can be decomposed into a set of time step problems, and then iterated between the time step problems to solve the decoupled problems until convergence.

Description:
RELATED APPLICATION 
     This is a Continuation-in Part application of U.S. patent application Ser. No. 14/185,024, “Method for Solving Quadratic Programs for Convex Sets with Linear Equalities by an Alternating Direction Method of Multipliers with Optimized Step Sizes,” filed by Raghunathan et al. on Feb. 20, 2014. 
    
    
     FIELD OF THE INVENTION 
     The invention relates generally to solving a multi-stage stochastic convex quadratic program (QP) for a convex set with linear equalities, and more particularly to solving the multi-stage stochastic convex QP by an Alternating Direction Method of Multipliers (ADMM). 
     BACKGROUND OF THE INVENTION 
     The Alternating Direction Method of Multipliers (ADMM), which is a variant of an augmented Lagrangian scheme, uses partial updates for dual variables to solve a quadratic program. The ADMM alternately solves for a first variable while holding a second variable fixed, and solving for the second variable while holding the first variable fixed, see Boyd et al. “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends in Machine Learning, vol. 3, no. 1, p. 1-122, 2011. Rather than iterating until convergence as in the Lagrangian-augmented problem, which is possibly complex, the ADMM performs iterations of alternating steps of updates on subsets of the variables. 
     The ADMM is frequently used for solving structured convex quadratic programs (QP) in applications, such as for example, compressive sensing, image processing, machine learning, distributed optimization, regularized estimation, semi-definite programming, and real-time control and signal processing applications, to name but a few. 
     Of particular interest to the invention is the application of ADMM to model predictive control (MPC) problems for application to systems with uncertain parameters. MPC is a method for controlling constrained dynamical systems, such as robots, automotive vehicles, spacecrafts, HVAC systems, processing machines, by repeatedly solving a finite horizon optimal control formulated from a mathematical model of the system dynamics, constraints, and a user specified cost function. Uncertainty arises due to the presence of disturbances that affect the system dynamics or due to mismatch between the model of the system dynamics and the true system dynamics. When probability information can be associated to the uncertain values, a stochastic system is obtained. 
     For linear models of the systems with a quadratic cost function and subject to linear constraints, the MPC finite horizon optimal control problem can be formulated as a parametric QP. At every control cycle, a specific QP is generated from the parametric quadratic program and the current values of the system state, and possibly the current reference. Then, the QP is solved, and the first part of the control input sequence is applied as control input. At the following control cycle, a new optimization problem is solved from the updated system state. MPC has been increasingly applied to systems with fast dynamics where the MPC is implemented in a low computational power embedded processors. 
     U.S. Pat. No. 8,600,525 discloses an active set algorithm that can be used for solving MPC generated QP problems. U.S. Pat. No. 7,328,074 discloses a method of providing an active-set algorithm wherein an initial guess for an optimization problem is provided from the solution of a previous optimization. U.S. Publication 20060282177 discloses an interior point algorithm for solving quadratic programs that arise in the context of model predictive control of gas turbine engines. 
     However, the computational effort per iteration of those methods can be prohibitive for application to large scale problems. The complexity of the operations that are performed, such as the solution of large scale linear systems, makes them infeasible for the type of computing hardware commonly used in real-time control and signal processing applications. 
     Methods such as gradient methods and accelerated gradient methods cannot easily handle linear inequality constraints. Low complexity fast optimization methods have been developed for MPC. U.S. Publication 20120281929 discloses a method for solving quadratic programs with non-negative constraints and a method to optimize such method for application to MPC. 
     A fast gradient algorithm is described for an application to MPC by Richter et al., “Computational Complexity Certification for Real-Time MPC With Input Constraints Based on the Fast Gradient Method,” IEEE Trans. Automat. Contr. 57(6): 1391-1403, 2012. 
     A Lagrangian method for MPC is described by Kögel et al., “Fast predictive control of linear, time-invariant systems using an algorithm based on the fast gradient method and augmented Lagrange multipliers,” CCA 2011: 780-785, 2011. 
     However, those methods are limited in the application by the types of constraints that can be handled, e.g., Richter et al., or can only handle input constraints, or need to perform division operations, e.g., in U.S. Publication 20120281929, which are time consuming in the computing hardware for MPC, or by complex iteration, e.g., Kögel et al. This is mainly due to the need of achieving the solution of the Lagrangian-augmented problem, which is complex due to the presence of constraints, before updating the Lagrange multipliers, multipliers at every update. 
     Thus, there is a need to provide a method that can solve large scale problems with small computational complexity per iteration, rapid convergence when problems are feasible, and quick detection of infeasibility. 
     SUMMARY OF THE INVENTION 
     The invention relates to finding solution of stochastic quadratic programs (StQPs) wherein constraints are linear and the objectives are quadratic and depend on uncertain parameters. 
     StQPs arise in the fields of financial planning, inventory management among others. Of particular interest is the application of quadratic stochastic programming to model predictive control of machines and manufactoring processes and plants. 
     The embodiments of the invention provide a method for solving a stochastic quadratic programs (StQP) for a convex set with a set of general linear equalities and inequalities by an alternating direction method of multipliers (ADMM). The invention relates to solution of Stochastic Quadratic Programs (StQPs) wherein the constraints are linear, objective is quadratic and the constraints and objective depend on uncertain parameters. 
     The method of the invention determines an optimal solution, or alternatively certifies that no solution to the problem exists. The embodiments also provide a method for optimizing a step size β for the ADMM, which achieves convergence with the least number of iterations. The embodiments also provide a method for accelerating the convergence of the ADMM using conjugate gradient (CG) method. 
     Of special interest is the application of ADMM to StQPs that are solved for Stochastic Model Predictive Control (StMPC), where inequalities represent constraints on states and controls of a dynamic system, and the equalities represent the equations of the system dynamics that couple the variables of the optimization problem. The stochasticity in MPC arises from uncertain parameters in the system dynamics arising from disturbances or plant-model mismatch when probability information on the uncertain parameters is available. A uniform probability on the uncertain values may be assumed to solve an MPC problem with uncertain parameters and no probability information by StMPC. 
     In StMPC, the uncertain parameters are typically sampled to obtain a finite set of realizations. The optimization is performed over the finite set of realization to obtain optimal control action that in expectation is optimal over all the realizations of the uncertain parameters. A StMPC problem formulation is an instance of a StQP. 
     Prior art methods are known for minimizing the QP while satisfying linear equalities and the convex set. For example, active set methods and interior point methods are the most common iterative methods. The embodiments use such methods to solve the StQP. 
     The embodiments of the invention overcome the difficulties of prior art methods by performing only simple operations in iterations, in contrast with U.S. Pat. Nos. 8,600,525, 7,328,074, an U.S. Publication 20060282177. The embodiments do not involve divisions as in U.S. Publication 20120281929, and perform low complexity iterations by alternatively updating subsets of Lagrange multipliers by single iterations, rather than updating the entire set of multipliers, thus operating only on simple constraints. 
     The embodiments provide the optimal selection for the step size β of the ADMM iterations, which minimizes the number of iterations performed by the method to a least number, and hence achieves the maximal rate of execution for a StMPC controller implemented in a given computing hardware. 
     In one embodiment, the StMPC problem is decomposed into two blocks. The first block corresponds to an equality constrained QP, and the second block corresponds to a projection onto the StMPC inequalities and anticipativity constraints. The StMPC problem can also be decomposed into a set of time step problems, and then iterated between the time step problems to solve the decoupled problems until convergence. 
     Furthermore, the embodiments allow for an early and effective detection of infeasibility, meaning that the method, on termination, certifies that no solution exists for the QP problem. The embodiments also described a solution that minimizes violations of the constraints. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  is a flow diagram of a method for solving a quadratic program by an alternating direction method of multipliers (ADMM) according to embodiments of the invention; 
         FIG. 2A  is a schematic of a convex set Y defined as a non-negative orthant that can be solved by embodiments of the invention; 
         FIG. 2B  is a schematic of a convex region Y defined as hyperplanes that can be solved by embodiments of the invention; 
         FIG. 3  is a block diagram of a controller and machine using the ADMM according to embodiments of the invention; 
         FIG. 4  is a block diagram of the controller using the ADMM according to embodiments of the invention; 
         FIGS. 5A and 5B  are schematics of scenario generation for stochastic MPC control according to embodiments of the invention; 
         FIG. 6  is a flow diagram of stochastic MPC based on scenario enumeration according to embodiments of the invention; and 
         FIG. 7  is a flow diagram for a CG process to accelerate convergence of ADMM according to embodiments of the invention. 
     
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
     Controller and Machine 
       FIG. 3  shows a controller  301  and a machine  302  using alternating direction method of multipliers (ADMM) according to embodiments of the invention. Of particular interest to the invention is the application of ADMM to model predictive control (MPC) problems with uncertain parameters, herein called the Stochastic MPC (StMPC). 
     The term “machine” is used generally because it is well understood that MPC has been used for decades in chemical and oil refineraries. Generally, models used in MPC are intended to represent the behavior of complex dynamical systems. The additional computations of the MPC is generally not needed to provide adequate control of simple systems, which are often controlled well by proportional integral-derivative PID controllers. 
     Common dynamic characteristics that are difficult for PID controllers include large time delays, constraints, multiple control inputs, and high-order dynamics. MPC models predict the change in the dependent variables of the modeled system that will be caused by changes in the independent variables. For example, in a chemical process, independent variables that can be adjusted by the controller are often either the setpoints of regulatory PID controllers (pressure, flow, temperature, etc.) or the final control element (valves, dampers, etc.). Independent variables that cannot be adjusted by the controller are used as disturbances. Dependent variables in these processes are other measurements that represent either control objectives or process constraints. 
     Stochastic Model Predictive Control (StMPC) 
     Stochastic model predictive control (StMPC) is a process for controlling systems subject to stochastic uncertainty. As shown for example in  FIG. 3 , the machine  302  that is to be controlled is subject to an uncertainty  307 , and receives an input  304 . The input is generated by the controller  301  so that the machine output  305  follows a desired behavior  303 , using feedback from a machine state  306  and an estimate of the uncertainty  309  obtained possibly by an uncertainty estimator  308 . In StMPC, the behavior of the machine subject to uncertainty is described by a linear stochastic difference Equation
 
 x ( t+ 1)= A   m ( w ( t )) x ( t )+ B   m ( w ( t )) u ( t )+ G   m ( w ( t ))
 
 p ( w ( t ))=ƒ p ( w ( t− 1), . . . , w ( t−t   w )),  (1)
 
where x∈             n     x    is the difference Equation state which contains the machine state and information of the desired behavior of the machine, u∈           n     u    is the input, w∈           n     w    is a stochastic vector that represents the uncertainty, possibly with a finite discrete set of values allowed, p is a probability distribution function, ƒ p  is a function that describes p, t w  is a length of time steps t that defines the probability distribution, and A m , B m , G m  are parameters of the model.

     The states and inputs of the system whose behavior is described in Equation (1) can be subject to constraints
 
 x ( t+k )∈ X   k|t   ,u ( t+k )∈ U   k|t ,  (2)
 
where X k|t , U k|t  are admissible sets for state and input, respectively. St MPC selects the control input by solving the optimal control problem formulated from Equations (1) and (2) as
 
                           ⁢         min     u   t       ⁢     E   ⁡     [         x     N   |   t     T     ⁢     Px     N   |   t         +       ∑     k   =   0     N     ⁢       x     k   |   t     T     ⁢     Qx     k   |   t           +       u     k   |   t     T     ⁢     Ru     k   |   t           ]         ⁢     
     ⁢       s   .   t   .           ⁢     x       k   +   1     |   t         =           A   m     ⁡     (     w   ⁡     (     t   +   k     )       )       ⁢     x     k   |   t         +         B   m     ⁡     (     w   ⁡     (     t   +   k     )       )       ⁢     u     k   |   t         +       G   m     ⁡     (     w   ⁡     (     t   +   k     )       )           ⁢     
     ⁢           ⁢       p   ⁡     (     w   ⁡     (     t   +   k     )       )       =       f   p     ⁡     (       w   ⁡     (     t   +   k   -   1     )       ,   •   ,     w   ⁡     (     t   +   k   -     t   w       )         )         ⁢     
     ⁢           ⁢         x     k   |   t       ∈     X     k   |   t         ,       u     k   |   t       ∈     U     k   |   t           ⁢     
     ⁢           ⁢         x     0   |   t       =     x   ⁡     (   t   )         ,               (   3   )               
where a k|t  is a predicted value of a generic vector a for k steps ahead of t, T is a transpose operator, Q and P are positive semidefinite matrices, R is a positive definite matrix, the positive integer N is a MPC prediction time horizon, and u t =(u 0|t , . . . , u N-1|t ) is the control policy.
 
     In the problem described by Equation (3), x is a random vector due to the effect of w. Hence, the problem is difficult to solve in such form. 
     Scenario-Enumeration Stochastic MPC Control 
     Scenario-enumeration stochastic MPC control operates as described for  FIG. 4 . The stochastic optimal control problem  406  is formulated from an uncertain machine model  405 , and a set of scenarios  402 , which are computed from uncertainty statistics  401 , and possibly from a current estimate of the uncertainty  309 . Scenario information  402  is used in conjunction with the cost function  403  to formulate the stochastic control problem  406 , which uses constraints  404 , the machine feedback  306  and the machine desired behavior  303 . Then, the stochastic optimal control problem  406  is solved during input computation  407 , from which the machine input  305  is determined. 
     Scenarios 
     The scenarios are sequences of possible realizations w, i.e.,
 
 s   k|t ( j )=[ w   0|t   j    . . . w   N-1|t   j ],
 
and to each scenario, a probability π(s k|t (j)) can be associated via p and ƒ p . The scenarios can be generated as described for  FIGS. 5A and 5B . From a current estimate w(t) of the uncertainty  501  r 1  possible values for the future uncertainty are selected w 0|t   j , j=1, . . . r 1  for instance based on their likelihood as described by the function ƒ p . These values are assigned to the nodes  502 . Then, the process is repeated from every node. A scenario is a sequence of connected nodes  503  and  504 . The scenarios can share an initial part, such as  505  shared between  503  and  504 . The scenarios can have different lengths, and the number of nodes generated from a previous node can change. A particular case is shown in  FIG. 5B  where r same length scenarios are generated from a single  506  node. The scenarios only share the initial node  506 .
 
     Stochastic MPC 
     The stochastic MPC based on scenario enumeration is described for  FIG. 6 . A current state  306  of the machine is read  601 , and a current estimate of the uncertainty is obtained  602 , the scenarios are generated  603 , the stochastic optimal control problem is constructed  604 , and solved to determined  605  the input. The computed input  304  is applied  606  to the machine and the cycle is repeated  607  when a new machine state is available. 
     In scenario-enumeration MPC for N r  scenarios, which is a subset of the possible scenarios on the MPC prediction horizon, the scenario-enumeration stochastic MPC problem is formulated as 
                       min       {     (       x   t   r     ,     u   t   r       )     }       r   =   1       N   r         ⁢       ∑     r   =   0       N   r       ⁢       π   ⁡     (       s     k   |   t       ⁡     (   r   )       )       ⁡     [           (     x     N   |   t     r     )     T     ⁢     Px     N   |   t     r       +       ∑     k   =   0       N   -   1       ⁢         (     x     k   |   t     r     )     T     ⁢     Qx     k   |   t     r         +         (     u     k   |   t     r     )     T     ⁢     Ru     k   |   t     r         ]           ⁢     
     ⁢           ⁢       s   .   t   .           ⁢     x       k   +   1     |   t         =           A   m     ⁡     (     w     k   |   t     r     )       ⁢     x     k   |   t     r       +         B   m     ⁡     (     w     k   |   t     r     )       ⁢     u     k   |   t     r       +       G   m     ⁡     (     w     k   |   t     r     )           ⁢     
     ⁢           ⁢         ∀   r     =   1     ,   …   ⁢           ,     N   r     ,     k   =   0     ,   …   ⁢           ,     N   -   1       ⁢     
     ⁢           ⁢         x       k   +   1     |   t     r     ∈     X       k   +   1     |   t         ,       u     k   |   t     r     ∈     U     k   |   t           ⁢     
     ⁢           ⁢         ∀   r     =   1     ,   …   ⁢           ,     N   r     ,     k   =   0     ,   …   ⁢           ,     N   -   1       ⁢     
     ⁢           ⁢         x     0   |   t     r     =         x   ⁡     (   t   )       ⁢     ∀   r       =   1       ,   …   ⁢           ,     N   r       ⁢     
     ⁢           ⁢         w     0   |   t     r     =         w   ⁡     (   t   )       ⁢     ∀   r       =   1       ,   …   ⁢           ,     N   r       ⁢     
     ⁢           ⁢         u     k   |   t         r   min     ⁡     (   k   )         =     u     k   |   t     r       ,     ∀     r   ≠       r   min     ⁡     (   k   )           ,       r   ∈       R     k   |   t       ⁢           ⁢   and   ⁢           ⁢       r   min     ⁡     (   k   )           =       min   ⁡     (     R     k   |   t       )       .                 (   4   )               
where x t   r =(x 1|t   r  . . . x N|t   r ), u t   r =(u 0|t   r  . . . u N-1|t   r ), and R k|t  is the set of all scenarios, which are equal from the beginning until the step k−1, that is,
 
 r,r′∈R   k|t               s   k′|t ( r )= s   k′|t ( r ′)∀ k′= 0, . . . , k− 1.

     Due to the scenario enumeration, all vectors in Equation (4) are now deterministic, so that Equation (4) can be formulated as 
                       min       {     y   r     }       r   =   1       N   r         ⁢       ∑     r   =   1       N   r       ⁢     (         1   2     ⁢     y   r   T     ⁢     Q   r     ⁢     y   r       +       q   r   T     ⁢     y   r         )         ⁢     
     ⁢           s   .   t   .           ⁢     A   r       ⁢     y   r       =         b   r     ⁢     ∀   r       =   1       ,   …   ⁢           ,     N   r       ⁢           ⁢           y   r     ∈       Y   r     ⁢     ∀   r         =   1     ,   …   ⁢           ,     N   r       ⁢     
     ⁢           y   _     k       r   min     ⁡     (   k   )         =         y   _     k   r     ⁢     ∀     r   ≠       r   min     ⁡     (   k   )               ,       r   ∈       R     k   |   t       ⁢           ⁢   and   ⁢           ⁢       r   min     ⁡     (   k   )           =     min   ⁡     (     R     k   |   t       )         ,             (   5   )               
where for each scenario r=1, . . . , N r , k max (r) is the largest index for which r∈R k     max     (r)|t . The variables  y   k   r  are called the non-anticipativity constrained variables. Note that, k max (r)≧0, which means that at a minimum the control at the current time t must be equated across all scenarios and
 
     
       
         
           
             
               
                 y 
                 r 
               
               = 
               
                 ( 
                 
                   
                     y 
                     r 
                   
                   , 
                   
                     
                       y 
                       _ 
                     
                     0 
                     r 
                   
                   , 
                   … 
                   ⁢ 
                   
                       
                   
                   , 
                   
                     
                       
                         y 
                         _ 
                       
                       k 
                       r 
                     
                     ⁢ 
                     
                       max 
                       
                         ( 
                         r 
                         ) 
                       
                     
                   
                 
                 ) 
               
             
             , 
             
               
                 y 
                 r 
               
               = 
               
                 ( 
                 
                   
                     x 
                     t 
                     r 
                   
                   = 
                   
                     u 
                     t 
                     r 
                   
                 
                 ) 
               
             
             , 
             
               
                 q 
                 r 
               
               = 
               0 
             
           
         
       
       
         
           
             
               Π 
               x 
               r 
             
             = 
             
               ( 
               
                 
                   
                     
                       π 
                       ⁡ 
                       
                         ( 
                         
                           
                             s 
                             
                               1 
                               | 
                               t 
                             
                           
                           ⁡ 
                           
                             ( 
                             r 
                             ) 
                           
                         
                         ) 
                       
                     
                   
                   
                     
                         
                     
                   
                   
                     0 
                   
                 
                 
                   
                     
                         
                     
                   
                   
                     ⋱ 
                   
                   
                     
                         
                     
                   
                 
                 
                   
                     0 
                   
                   
                     … 
                   
                   
                     
                       π 
                       ⁡ 
                       
                         ( 
                         
                           
                             s 
                             
                               
                                 N 
                                 - 
                                 1 
                               
                               | 
                               t 
                             
                           
                           ⁡ 
                           
                             ( 
                             r 
                             ) 
                           
                         
                         ) 
                       
                     
                   
                 
               
               ) 
             
           
         
       
       
         
           
             
               Π 
               u 
               r 
             
             = 
             
               ( 
               
                 
                   
                     
                       π 
                       ⁡ 
                       
                         ( 
                         
                           
                             s 
                             
                               0 
                               | 
                               t 
                             
                           
                           ⁡ 
                           
                             ( 
                             r 
                             ) 
                           
                         
                         ) 
                       
                     
                   
                   
                     
                         
                     
                   
                   
                     0 
                   
                 
                 
                   
                     
                         
                     
                   
                   
                     ⋱ 
                   
                   
                     
                         
                     
                   
                 
                 
                   
                     0 
                   
                   
                     … 
                   
                   
                     
                       π 
                       ⁡ 
                       
                         ( 
                         
                           
                             s 
                             
                               
                                 N 
                                 - 
                                 1 
                               
                               | 
                               t 
                             
                           
                           ⁡ 
                           
                             ( 
                             r 
                             ) 
                           
                         
                         ) 
                       
                     
                   
                 
               
               ) 
             
           
         
       
       
         
           
             
               Q 
               r 
             
             = 
             
               ( 
               
                 
                   
                     
                       
                         Π 
                         x 
                         r 
                       
                       ⊗ 
                       Q 
                     
                   
                   
                     0 
                   
                   
                     0 
                   
                   
                     0 
                   
                 
                 
                   
                     0 
                   
                   
                     
                       
                         π 
                         ⁡ 
                         
                           ( 
                           
                             
                               s 
                               
                                 k 
                                 | 
                                 t 
                               
                             
                             ⁡ 
                             
                               ( 
                               r 
                               ) 
                             
                           
                           ) 
                         
                       
                       ⁢ 
                       P 
                     
                   
                   
                     0 
                   
                   
                     0 
                   
                 
                 
                   
                     0 
                   
                   
                     0 
                   
                   
                     
                       
                         Π 
                         u 
                         r 
                       
                       ⊗ 
                       R 
                     
                   
                   
                     0 
                   
                 
                 
                   
                     0 
                   
                   
                     0 
                   
                   
                     0 
                   
                   
                     0 
                   
                 
               
               ) 
             
           
         
       
       
         
           
             b 
             = 
             
               ( 
               
                 
                   
                     
                       
                         
                           
                             A 
                             m 
                           
                           ⁡ 
                           
                             ( 
                             
                               w 
                               
                                 0 
                                 | 
                                 t 
                               
                               r 
                             
                             ) 
                           
                         
                         ⁢ 
                         
                           x 
                           ⁡ 
                           
                             ( 
                             t 
                             ) 
                           
                         
                       
                       + 
                       
                         
                           G 
                           m 
                         
                         ⁡ 
                         
                           ( 
                           
                             w 
                             
                               0 
                               | 
                               t 
                             
                             r 
                           
                           ) 
                         
                       
                     
                   
                   
                     
                         
                     
                   
                 
                 
                   
                     
                       
                         G 
                         m 
                       
                       ⁡ 
                       
                         ( 
                         
                           w 
                           
                             1 
                             | 
                             t 
                           
                           r 
                         
                         ) 
                       
                     
                   
                   
                     
                         
                     
                   
                 
                 
                   
                     ⋮ 
                   
                   
                     
                         
                     
                   
                 
                 
                   
                     
                       
                         G 
                         m 
                       
                       ⁡ 
                       
                         ( 
                         
                           w 
                           
                             
                               N 
                               - 
                               1 
                             
                             | 
                             t 
                           
                           r 
                         
                         ) 
                       
                     
                   
                   
                     
                         
                     
                   
                 
                 
                   
                     
                       
                         
                           
                             0 
                           
                         
                         
                           
                             ⋮ 
                           
                         
                         
                           
                             0 
                           
                         
                       
                       } 
                     
                   
                   
                     
                       
                         
                           n 
                           u 
                         
                         ⁡ 
                         
                           ( 
                           
                             
                               
                                 k 
                                 max 
                               
                               ⁡ 
                               
                                 ( 
                                 r 
                                 ) 
                               
                             
                             + 
                             1 
                           
                           ) 
                         
                       
                       ⁢ 
                       times 
                     
                   
                 
               
               ) 
             
           
         
       
       
         
           
             
               Y 
               r 
             
             = 
             
               
                 
                   
                     X 
                     
                       1 
                       | 
                       t 
                     
                   
                   × 
                   … 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     X 
                     
                       N 
                       | 
                       t 
                     
                   
                   × 
                   
                     U 
                     
                       0 
                       | 
                       t 
                     
                   
                   × 
                   … 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     U 
                     
                       
                         N 
                         - 
                         1 
                       
                       | 
                       t 
                     
                   
                 
                 
                   ︸ 
                   
                     = 
                     
                       : 
                       
                         Y 
                         r 
                       
                     
                   
                 
               
               × 
               
                 
                   R 
                   
                     
                       n 
                       u 
                     
                     ⁡ 
                     
                       ( 
                       
                         
                           
                             k 
                             max 
                           
                           ⁡ 
                           
                             ( 
                             r 
                             ) 
                           
                         
                         + 
                         1 
                       
                       ) 
                     
                   
                 
                 . 
               
             
           
         
       
     
     The matrix A r  is defined so that, 
     
       
         
           
             
                 
             
             ⁢ 
             
               
                 A 
                 r 
               
               = 
               
                 ( 
                 
                   
                     
                       
                         
                           A 
                           ^ 
                         
                         r 
                       
                     
                     
                       
                         
                           B 
                           ^ 
                         
                         r 
                       
                     
                     
                       0 
                     
                   
                   
                     
                       0 
                     
                     
                       
                         ( 
                         
                           
                             
                               
                                 - 
                                 
                                   I 
                                   
                                     
                                       n 
                                       u 
                                     
                                     ⁡ 
                                     
                                       ( 
                                       
                                         
                                           
                                             k 
                                             max 
                                           
                                           ⁡ 
                                           
                                             ( 
                                             r 
                                             ) 
                                           
                                         
                                         + 
                                         1 
                                       
                                       ) 
                                     
                                   
                                 
                               
                             
                             
                               0 
                             
                           
                         
                         ) 
                       
                     
                     
                       
                         I 
                         
                           
                             n 
                             u 
                           
                           ⁡ 
                           
                             ( 
                             
                               
                                 
                                   k 
                                   max 
                                 
                                 ⁡ 
                                 
                                   ( 
                                   r 
                                   ) 
                                 
                               
                               + 
                               1 
                             
                             ) 
                           
                         
                       
                     
                   
                 
                 ) 
               
             
           
         
       
       
         
           
             
                 
             
             ⁢ 
             
               where 
               , 
               
                 
 
               
               ⁢ 
               
                   
               
               ⁢ 
               
                 
                   
                     A 
                     ^ 
                   
                   r 
                 
                 = 
                 
                   ( 
                   
                     
                       
                         
                           I 
                           
                             n 
                             x 
                           
                         
                       
                       
                         0 
                       
                       
                         … 
                       
                       
                         0 
                       
                       
                         0 
                       
                     
                     
                       
                         
                           - 
                           
                             
                               A 
                               m 
                             
                             ⁡ 
                             
                               ( 
                               
                                 w 
                                 
                                   1 
                                   | 
                                   t 
                                 
                                 r 
                               
                               ) 
                             
                           
                         
                       
                       
                         
                           I 
                           
                             n 
                             x 
                           
                         
                       
                       
                         … 
                       
                       
                         0 
                       
                       
                         0 
                       
                     
                     
                       
                         ⋮ 
                       
                       
                         ⋮ 
                       
                       
                         ⋱ 
                       
                       
                         ⋮ 
                       
                       
                         ⋮ 
                       
                     
                     
                       
                         0 
                       
                       
                         0 
                       
                       
                         … 
                       
                       
                         
                           I 
                           
                             n 
                             x 
                           
                         
                       
                       
                         0 
                       
                     
                     
                       
                         0 
                       
                       
                         0 
                       
                       
                         … 
                       
                       
                         
                           - 
                           
                             
                               A 
                               m 
                             
                             ⁡ 
                             
                               ( 
                               
                                 w 
                                 
                                   
                                     N 
                                     - 
                                     1 
                                   
                                   | 
                                   t 
                                 
                                 r 
                               
                               ) 
                             
                           
                         
                       
                       
                         
                           I 
                           
                             n 
                             x 
                           
                         
                       
                     
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
                 
             
             ⁢ 
             
               and 
               , 
               
                 
 
               
               ⁢ 
               
                   
               
               ⁢ 
               
                 
                   
                     B 
                     ^ 
                   
                   r 
                 
                 = 
                 
                   ( 
                   
                     
                       
                         
                           - 
                           
                             
                               B 
                               m 
                             
                             ⁡ 
                             
                               ( 
                               
                                 w 
                                 
                                   0 
                                   | 
                                   t 
                                 
                                 r 
                               
                               ) 
                             
                           
                         
                       
                       
                         … 
                       
                       
                         0 
                       
                     
                     
                       
                         
                             
                         
                       
                       
                         ⋱ 
                       
                       
                         
                             
                         
                       
                     
                     
                       
                         0 
                       
                       
                         … 
                       
                       
                         
                           - 
                           
                             
                               B 
                               m 
                             
                             ⁡ 
                             
                               ( 
                               
                                 w 
                                 
                                   
                                     N 
                                     - 
                                     1 
                                   
                                   | 
                                   t 
                                 
                                 r 
                               
                               ) 
                             
                           
                         
                       
                     
                   
                   ) 
                 
               
               , 
               
                 
 
               
               ⁢ 
               
                 
                   I 
                   
                     
                       n 
                       u 
                     
                     ⁡ 
                     
                       ( 
                       
                         
                           
                             k 
                             max 
                           
                           ⁡ 
                           
                             ( 
                             r 
                             ) 
                           
                         
                         + 
                         1 
                       
                       ) 
                     
                   
                 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 is 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 an 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 identity 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 matrix 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 of 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 size 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 
                   
                     n 
                     u 
                   
                   ⁡ 
                   
                     ( 
                     
                       
                         
                           k 
                           max 
                         
                         ⁡ 
                         
                           ( 
                           r 
                           ) 
                         
                       
                       + 
                       1 
                     
                     ) 
                   
                 
                 × 
                 
                   
                     
                       n 
                       u 
                     
                     ⁡ 
                     
                       ( 
                       
                         
                           
                             k 
                             max 
                           
                           ⁡ 
                           
                             ( 
                             r 
                             ) 
                           
                         
                         + 
                         1 
                       
                       ) 
                     
                   
                   . 
                 
               
             
           
         
       
     
     Scenario Decomposition-Based ADMM 
     As shown in  FIG. 1 , the embodiments of the invention provide a method for solving a quadratic program (QP) for a convex set with general linear equalities and inequalities by an alternating direction method of multipliers (ADMM). The method optimizes the step size β of the ADMM to minimize the number of required iterations thus minimizing the amount of time of the processor needed to obtain the solution. 
     In one embodiment of the invention, the quadratic program (5) resulting from scenario-enumeration stochastic MPC problem can be reformulated as, 
                       min       {     y   r     }       r   =   1       N   r         ⁢       ∑     r   =   1       N   r       ⁢     (         1   2     ⁢     y   r   T     ⁢     Q   r     ⁢     y   r       +       q   r   T     ⁢     y   r         )         ⁢     
     ⁢           s   .   t   .           ⁢     A   r       ⁢     y   r       =         b   r     ⁢     ∀   r       =   1       ,   …   ⁢           ,     N   r     ,     
     ⁢     y   ∈     Y   sd                 (   6   )               
where a convex set Y sd  is defined as
 
     
       
         
           
             
               
                 
                   
                     Y 
                     sd 
                   
                   = 
                   
                     
                       { 
                       
                         
                           ( 
                           
                             
                               y 
                               1 
                             
                             , 
                             … 
                             ⁢ 
                             
                                 
                             
                             , 
                             
                               y 
                               
                                 N 
                                 r 
                               
                             
                           
                           ) 
                         
                         | 
                         
                           
                             
                               
                                 
                                   
                                     
                                       y 
                                       r 
                                     
                                     ∈ 
                                     
                                       
                                         Y 
                                         r 
                                       
                                       ⁢ 
                                       
                                         ∀ 
                                         r 
                                       
                                     
                                   
                                   = 
                                   1 
                                 
                                 , 
                                 … 
                                 ⁢ 
                                 
                                     
                                 
                                 , 
                                 
                                   N 
                                   r 
                                 
                               
                             
                           
                           
                             
                               
                                 
                                   
                                     
                                       
                                         
                                           
                                             
                                               y 
                                               _ 
                                             
                                             k 
                                             
                                               r 
                                               min 
                                             
                                           
                                           ⁡ 
                                           
                                             ( 
                                             k 
                                             ) 
                                           
                                         
                                         = 
                                         
                                           
                                             
                                               y 
                                               _ 
                                             
                                             
                                               
                                                 r 
                                                 min 
                                               
                                               ⁡ 
                                               
                                                 ( 
                                                 k 
                                                 ) 
                                               
                                             
                                           
                                           ⁢ 
                                           
                                             ∀ 
                                             
                                               r 
                                               ≠ 
                                               
                                                 r 
                                                 k 
                                               
                                             
                                           
                                         
                                       
                                       , 
                                       
                                         r 
                                         ∈ 
                                         
                                           R 
                                           
                                             k 
                                             | 
                                             t 
                                           
                                         
                                       
                                     
                                   
                                 
                                 
                                   
                                     
                                       
                                         and 
                                         ⁢ 
                                         
                                             
                                         
                                         ⁢ 
                                         
                                           
                                             r 
                                             min 
                                           
                                           ⁡ 
                                           
                                             ( 
                                             k 
                                             ) 
                                           
                                         
                                       
                                       = 
                                       
                                         min 
                                         ⁡ 
                                         
                                           ( 
                                           
                                             R 
                                             
                                               k 
                                               | 
                                               t 
                                             
                                           
                                           ) 
                                         
                                       
                                     
                                   
                                 
                               
                             
                           
                         
                       
                       } 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   7 
                   ) 
                 
               
             
           
         
       
     
     Convex Sets 
       FIGS. 2A and 2B  show example convex sets Y that can be solved by the embodiments of the invention. In a Euclidean space, a set is convex when, for every pair of points  201 - 202  in the set, every point on a straight line segment  210  that joins the pair of points is also in the set. Specifically,  FIG. 2A  shows a convex set Y  200  defined as a non-negative orthant  200 .  FIG. 2B  is a schematic of a convex set Y defined as hyperplanes  220 . Observe that non-anticipativity constraints are included in the convex set Y sd . 
     Consider reformulating Equation (6) as, 
                       min     y   ,   w       ⁢       ∑     r   =   1       N   r       ⁢     (         1   2     ⁢     y   r   T     ⁢     Q   r     ⁢     y   r       +       q   r   T     ⁢     y   r         )         ⁢     
     ⁢         s   .   t   .           ⁢     A   r       ⁢     y   r       =         b   r     ⁢     ∀   r       =   1       ,   …   ⁢           ,     N   r     ,     
     ⁢     y   =   w     ,     w   ∈     Y   sd               (   8   )               
where y=(y 1 , . . . , y N     r   )∈             n  with n=N r (n x +n u +n u (k max (r)+1)). The variables w∈           n  are required to be in the convex set Y SD . The advantage of Equation (8) is that the inequalities and the scenario coupling non-anticipativity constraints are placed on separate variables w, coupled with the others by y=w. The variables y is the linear subspace constrained variables, and w will be called the convex set constrained variables.

     The ADMM procedure dualizes the constraints y r =w r  into the objective function using scaled multipliers β r λ r  respectively where β r &gt;0. Additionally, the objective is also augmented with a penalty on the squared norm of the violation of the dualized equality constraints. Thus, we obtain 
                       min     y   ,   w       ⁢       ∑     r   =   1       N   r       ⁢       L   r     ⁡     (       y   r     ,     w   r     ,     λ   r       )           ⁢     
     ⁢         s   .   t   .           ⁢     A   r       ⁢     y   r       =         b   r     ⁢     ∀   r       =   1       ,   …   ⁢           ,     N   r     ,     
     ⁢     w   ∈     Y   sd               (   9   )               
where
 
     
       
         
           
             
               
                 
                   
                     
                       L 
                       r 
                     
                     ⁡ 
                     
                       ( 
                       
                         
                           y 
                           r 
                         
                         , 
                         
                           w 
                           r 
                         
                         , 
                         
                           λ 
                           r 
                         
                       
                       ) 
                     
                   
                   := 
                   
                     
                       
                         1 
                         2 
                       
                       ⁢ 
                       
                         y 
                         r 
                         T 
                       
                       ⁢ 
                       
                         Q 
                         r 
                       
                       ⁢ 
                       
                         y 
                         r 
                       
                     
                     + 
                     
                       
                         q 
                         r 
                         T 
                       
                       ⁢ 
                       
                         y 
                         r 
                       
                     
                     + 
                     
                       
                         
                           β 
                           r 
                         
                         2 
                       
                       ⁢ 
                       
                         
                           
                              
                             
                               
                                 y 
                                 r 
                               
                               - 
                               
                                 w 
                                 r 
                               
                               - 
                               
                                 λ 
                                 r 
                               
                             
                              
                           
                           2 
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   10 
                   ) 
                 
               
             
           
         
       
     
     In Equation (9), w and y are coupled only by the objective function. The steps  110 ,  120 ,  130  and  140  in the ADMM procedure respectively are: 
     
       
         
           
             
               
                 
                   
                     
                       y 
                       
                         l 
                         + 
                         1 
                       
                     
                     = 
                     
                       arg 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         
                           min 
                           y 
                         
                         ⁢ 
                         
                           
                             ∑ 
                             
                               r 
                               = 
                               1 
                             
                             
                               N 
                               r 
                             
                           
                           ⁢ 
                           
                             
                               L 
                               r 
                             
                             ⁡ 
                             
                               ( 
                               
                                 
                                   y 
                                   r 
                                 
                                 , 
                                 
                                   w 
                                   r 
                                   l 
                                 
                                 , 
                                 
                                   λ 
                                   r 
                                   l 
                                 
                               
                               ) 
                             
                           
                         
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       
                         
                           s 
                           . 
                           t 
                           . 
                           
                               
                           
                           ⁢ 
                           
                             A 
                             r 
                           
                         
                         ⁢ 
                         
                           y 
                           r 
                         
                       
                       = 
                       
                         
                           
                             b 
                             r 
                           
                           ⁢ 
                           
                             ∀ 
                             r 
                           
                         
                         = 
                         1 
                       
                     
                     , 
                     … 
                     ⁢ 
                     
                         
                     
                     , 
                     
                       N 
                       r 
                     
                   
                 
               
               
                 
                   ( 
                   11 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       w 
                       
                         l 
                         + 
                         1 
                       
                     
                     = 
                     
                       arg 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         
                           min 
                           w 
                         
                         ⁢ 
                         
                           
                             ∑ 
                             
                               r 
                               = 
                               1 
                             
                             
                               N 
                               r 
                             
                           
                           ⁢ 
                           
                             
                               L 
                               r 
                             
                             ⁡ 
                             
                               ( 
                               
                                 
                                   y 
                                   r 
                                   
                                     l 
                                     + 
                                     1 
                                   
                                 
                                 , 
                                 w 
                                 , 
                                 
                                   
                                     λ 
                                     ~ 
                                   
                                   r 
                                   l 
                                 
                               
                               ) 
                             
                           
                         
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       s 
                       . 
                       t 
                       . 
                       
                           
                       
                       ⁢ 
                       w 
                     
                     ∈ 
                     
                       Y 
                       sd 
                     
                   
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
             
               
                 
                   
                     λ 
                     
                       l 
                       + 
                       1 
                     
                   
                   = 
                   
                     
                       λ 
                       l 
                     
                     + 
                     
                       w 
                       
                         l 
                         + 
                         1 
                       
                     
                     - 
                     
                       
                         y 
                         
                           l 
                           + 
                           1 
                         
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
       
     
     The update step  140  of λ in Equation (13) scales linearly in the number of scenarios. We show in the following that Equations (11) and (12) also decouple by scenarios. 
     Observe that the objective functions and constraints in Equation (11) are decoupled by scenarios. Hence, the update  111  can be rewritten as,
 
 y   r   l+1   =M   r ( w   r   l +λ r   l   −q   r /β r )+ N   r   b   r  
 
where,
 
 M   r   :=Z   r ( Z   r   T ( Q   r /β r   +I   n ) Z   r ) −1   Z   r   T ,
 
 N   r :=( I   n   −M   r   Q   r /β r ) R   r ( A   r   R   r ) −1 ,  (14)
 
with R r ,Z r  denote an orghonormal basis for the range space of A r   T  and null space of A r , respectively. Thus, the update for y decouples by scenario and scales linearly with the number of scenarios.
 
     In (12), the objective function is component wise separable in the w. Using w r =(w r , w   r ) the constraint set w∈Y sd  can be posed as,
 
 w   r   ∈Y   r   ∀r= 1, . . . , n   r  and
 
   w     k   r     min       (k)     = w     k   r   ∀r≠r   min ( k ), r∈R   k|t  and  r   min ( k )=min( R   k|t )  (15)
 
     Note that  w   r  are only constrained by the equality of the values across scenarios and are not limited by bounds. Hence, the update  121  in (12) can be obtained as 
                       w   r     l   +   1       =       P   Y     ⁡     (         (     y   r     )       l   +   1       -       (     λ   r     )     l       )         ⁢     
     ⁢         (       w   _     k   r     )       l   +   1       =         1          R     k   |   t              ⁢       ∑       r   ′     ∈     R     k   |   t           ⁢       (         (       y   _     k     r   ′       )       l   +   1       -       (       λ   _     r     )     l       )     ⁢     ∀   k           =       ,   0   ,   …   ⁢           ,     N   -   1     ,     r   ∈     R     k   |   t                 (   16   )               
where we have used λ r =(λ r , λ   r ) and P Y (x) denotes the projection of vector x onto the set Y. Even though the set Y sd  is not simple we can still compute the update for w in a manner that scales linearly in the number of scenarios.
 
     Calling ADMM-CG 
     A check is made  133  to call the CG process every n admm  iterations, where n admm  is predetermined. If yes, then the CG process  700  is called  135 . The ADMM-CG is explained in detail later in this invention. Otherwise, it is proceeded to step  150 . The ADMM-CG procedure is described in greater detail later in this invention with reference to  FIG. 7   
     Termination Condition 
     Step  150  checks if an optimality termination holds. Given a predetermined tolerance ∈, a termination condition  151  for an optimal feasible solution  160  is
 
 r   opt =max(∥ w   l+1   −w   l ∥,∥λ l+1 −λ l ∥)&lt;∈.  (17)
 
     If the termination condition for the feasible solution is satisfied, then the optimal feasible solution is output  160 . The termination condition  151  in Equation (17) checks for the satisfaction to the tolerance greater ∈ greater than zero of a maximum of a norm of a change in the set constrained variable w from a current value to a value at a previous iteration, and a norm of the change in the Lagrange multiplier λ from a current value to a value at a previous iteration. 
     Otherwise, for ∈&gt;0, a termination condition  171  for certifying a solution to the problem is infeasibility is checked in step  170   
                     r   inf     =       max   ⁡     (              y     l   +   1       -     y   l            ,            w     l   +   1       -     w   l            ,     1   -           (       w     l   +   1       -     y     l   +   1         )     T     ⁢     λ     l   +   1                    w     l   +   1       -     y     l   +   1              ⁢          λ     l   +   1                    )       &lt;       ɛ   ⁢     
     ⁢           (       y     l   +   1       -     w     l   +   1         )     ∘     λ     l   +   1         ≥   0             (   18   )               
where ∘ denotes the element-wise multiplication of two vectors. If the termination condition for the infeasible solution is satisfied, then certification that the solution is infeasible can be signaled  180 .
 
     The termination condition  171  in Equation (18) is checked for the satisfaction of four conditions. 
     The first condition is a satisfaction to a tolerance ∈ greater than zero of a maximum of a normed change in the set constrained variable vector w from a current value to a value at a previous iteration. 
     The second condition is the satisfaction to a tolerance ∈ greater than zero of the normed change in the linear subspace constrained variable vector y from the current value to the value at the previous iteration. 
     The third condition checks for a deviation from 0 to a tolerance of not more than ∈ of an angle between the Lagrange multiplier vector λ and the vector resulting from a difference of the linear subspace constrained variable vector and the set constrained variable vector, i.e., (y−w), at the current value. 
     The fourth condition requires that a difference of the linear subspace constrained variable vector and the set constrained variable vector at the current value have, element-wise, an identical sign as the Lagrange multiplier vector. Otherwise, update  140  l=l+1, and perform the next iteration. 
     The general method can be implemented in a processor or other hardware as describe above connected to memory and input/output interfaces by buses as known in the art. 
     Optimal Parameter Choice 
     The choice of the step size β r , which ensures that a least number of iterations are required for termination of the ADMM method is
 
β r   opt =√{square root over (λ min ( Z   r   T   Q   r   Z   r )λ max ( Z   r   T   Q   r   Z   r ))},  (19)
 
where λ min (•),λ max (•) are minimal and maximal eigenvalues of contained matrix arguments. In other words, the optimal step size is determined as a square root of a product of minimum and maximum eigenvalues of a Hessian matrix of the scenarios problem pre- and post multiplied by an orthonormal basis for a null space of a linear equality constraint matrix.
 
     In another embodiment of the invention, a single value of β is chosen for all scenarios β r =β. In this case, the optimal value of the parameter is prescribed as,
 
β opt =√{square root over (λ min ( Z   T   QZ )λ max ( Z   T   QZ ))},  (20)
 
where, Z is the orthonormal basis for the vectors satisfying the constraints,
 
 A   r   y   r =0∀ r= 1, . . . , N   r  
 
   y     k   r     min       (k)     = y     k   r   ∀r≠r   k   ,r∈R   k|t  and  r   min ( k )=min( R   k|t )
 
and,
 
     
       
         
           
             Q 
             = 
             
               
                 ( 
                 
                   
                     
                       
                         Q 
                         1 
                       
                     
                     
                       
                           
                       
                     
                     
                       0 
                     
                   
                   
                     
                       
                           
                       
                     
                     
                       ⋱ 
                     
                     
                       
                           
                       
                     
                   
                   
                     
                       0 
                     
                     
                       … 
                     
                     
                       
                         Q 
                         
                           N 
                           r 
                         
                       
                     
                   
                 
                 ) 
               
               . 
             
           
         
       
     
     Scenario Decomposition and Conjugate Gradient (CG)-Based ADMM 
       FIG. 7  describes the steps involved in the ADMM-CG process  700  used to accelerate the convergence of the ADMM iterations. For the purposes of simplifying this description, the following are assumed. It is understood that the techniques described herein can be extended identically to multi-stage stochastic programs: Y r =[ y   r , y   r ], where  y   r  is a vector of lower bounds and  y   r  is vector of upper bounds for y r , and the scenarios are such that R k|t =Ø∀k&gt;1. In other words, the MPC problem is a two-stage stochastic program. 
     Hence,
 
 y   r =( y   r   , y     0   r ), w   r =( w   r   , w     0   r ),λ r =(λ r , λ   0   r )∉ r= 1, . . . , N   r ,
 
where
 
( y   +   ,w   + ,λ + )=ADMM( y,w ,λ)
 
denotes one iteration of ADMM, that is application of Equations (11-13), where (y k , w k , λ k )=(y, w, λ), and the parameter β r =β. For convenience, let v l =y l −λ l−1 .
 
     After every n admm  iterations of the ADMM-CG process  700 , as checked in step  133 , defined by steps in Equation (11-13), the ADMM-CG process is called  135 , see  FIG. 1 . 
     Identifying Linear Inequalities Expected to Hold as Equalities at Solution 
     Given an ADMM iterate (y l , w l , λ l )  705 , we define the index sets  710  as 
                           I   _     r     =     {       i   |       (     w   i   r     )     l       =           y   _     i   r     ⁢           ⁢   or   ⁢           ⁢       (     λ   i   r     )     l       ≥     ε   act         }       ⁢     
     ⁢       I   _     r     =     {       i   |       (     w   i   r     )     l       =           y   _     i   r     ⁢           ⁢   or   ⁢           ⁢       (     λ   i   r     )     l       ≤     -     ε   act           }       ⁢     
     ⁢         y   ^       r   ,   i       =     {                   y   _     i   r     ⁢           ⁢   if   ⁢           ⁢   i     ∈       I   _     r                       y   _     i   r     ⁢           ⁢   if   ⁢           ⁢   i     ∈       I   _     r                 0   ⁢           ⁢   otherwise           ,                 (   21   )               
where ∈ act  is a tolerance for estimating the inequality constraints in the convex set that are expected to hold as equalities at the solution to the StQP, and E r  is a matrix that extracts the components of w r  that are in  I   r ∪Ī r , and E r0  is a matrix that extracts the components of w r  that correspond to  y   0   r . The above active index sets in  I   r ,Ī r  are used as estimates of the set of indices of y r  that are on a bound for an optimal solution to Equation (6). The CG method  715  is used to solve the following system of linear equations called a CG-linear system:
 
                       M   ~     ⁡     (           λ   ^               λ   _           )       =     b   ~             (   22   )               
where
 
     
       
         
           
             
               
                 
                   
                     
                       
                         λ 
                         ^ 
                       
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 
                                   λ 
                                   ^ 
                                 
                                 1 
                               
                             
                           
                           
                             
                               ⋮ 
                             
                           
                           
                             
                               
                                 
                                   λ 
                                   ^ 
                                 
                                 
                                   N 
                                   r 
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                     , 
                     
                       
                         λ 
                         ^ 
                       
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 
                                   λ 
                                   ~ 
                                 
                                 2 
                               
                             
                           
                           
                             
                               ⋮ 
                             
                           
                           
                             
                               
                                 
                                   λ 
                                   ~ 
                                 
                                 
                                   N 
                                   r 
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                     , 
                     
                       q 
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 q 
                                 1 
                               
                             
                           
                           
                             
                               ⋮ 
                             
                           
                           
                             
                               
                                 q 
                                 
                                   N 
                                   r 
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                     , 
                     
                       b 
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 b 
                                 1 
                               
                             
                           
                           
                             
                               ⋮ 
                             
                           
                           
                             
                               
                                 b 
                                 
                                   N 
                                   r 
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                     , 
                     
                       
 
                     
                     ⁢ 
                     
                       Q 
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 Q 
                                 1 
                               
                             
                             
                               
                                   
                               
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               
                                   
                               
                             
                             
                               ⋱ 
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               
                                   
                               
                             
                             
                               
                                   
                               
                             
                             
                               
                                 Q 
                                 
                                   N 
                                   r 
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                     , 
                     
                       A 
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 A 
                                 1 
                               
                             
                             
                               
                                   
                               
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               
                                   
                               
                             
                             
                               ⋱ 
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               
                                   
                               
                             
                             
                               
                                   
                               
                             
                             
                               
                                 A 
                                 
                                   N 
                                   r 
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                     , 
                     
                       
 
                     
                     ⁢ 
                     
                       E 
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 E 
                                 1 
                               
                             
                             
                               
                                   
                               
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               
                                   
                               
                             
                             
                               ⋱ 
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               
                                   
                               
                             
                             
                               
                                   
                               
                             
                             
                               
                                 E 
                                 
                                   N 
                                   r 
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                     , 
                     
                       
                         E 
                         0 
                       
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 E 
                                 10 
                               
                             
                             
                               
                                 - 
                                 
                                   E 
                                   20 
                                 
                               
                             
                             
                               0 
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               
                                 E 
                                 10 
                               
                             
                             
                               
                                 
                                     
                                 
                                 ⁢ 
                                 0 
                               
                             
                             
                               
                                 - 
                                 
                                   E 
                                   30 
                                 
                               
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               ⋮ 
                             
                             
                               
                                   
                               
                             
                             
                               
                                   
                               
                             
                             
                               ⋱ 
                             
                           
                           
                             
                               
                                 E 
                                 10 
                               
                             
                             
                               
                                   
                               
                             
                             
                               0 
                             
                             
                               
                                 - 
                                 
                                   E 
                                   
                                     
                                       n 
                                       r 
                                     
                                     ⁢ 
                                     0 
                                   
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       Z 
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 Z 
                                 1 
                               
                             
                             
                               
                                   
                               
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               
                                   
                               
                             
                             
                               ⋱ 
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               
                                   
                               
                             
                             
                               
                                   
                               
                             
                             
                               
                                 Z 
                                 
                                   N 
                                   r 
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                     , 
                     
                       R 
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 R 
                                 1 
                               
                             
                             
                               
                                   
                               
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               
                                   
                               
                             
                             
                               ⋱ 
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               
                                   
                               
                             
                             
                               
                                   
                               
                             
                             
                               
                                 R 
                                 
                                   N 
                                   r 
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                     , 
                     
                       
 
                     
                     ⁢ 
                     
                       
                         Q 
                         ~ 
                       
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 
                                   ( 
                                   
                                     
                                       Z 
                                       1 
                                       T 
                                     
                                     ⁢ 
                                     
                                       Q 
                                       r 
                                     
                                     ⁢ 
                                     
                                       Z 
                                       1 
                                     
                                   
                                   ) 
                                 
                                 
                                   - 
                                   1 
                                 
                               
                             
                             
                               
                                   
                               
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               
                                   
                               
                             
                             
                               
                                 
                                     
                                 
                                 ⁢ 
                                 ⋱ 
                               
                             
                             
                               
                                   
                               
                             
                           
                           
                             
                               
                                   
                               
                             
                             
                               
                                   
                               
                             
                             
                               
                                 
                                   ( 
                                   
                                     
                                       Z 
                                       
                                         N 
                                         r 
                                       
                                       T 
                                     
                                     ⁢ 
                                     
                                       Q 
                                       
                                         N 
                                         r 
                                       
                                     
                                     ⁢ 
                                     
                                       Z 
                                       
                                         N 
                                         r 
                                       
                                     
                                   
                                   ) 
                                 
                                 
                                   - 
                                   1 
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       
                         M 
                         ~ 
                       
                       = 
                       
                         
                           ( 
                           
                             
                               
                                 E 
                               
                             
                             
                               
                                 
                                   E 
                                   0 
                                 
                               
                             
                           
                           ) 
                         
                         ⁢ 
                         Z 
                         ⁢ 
                         
                           Q 
                           ~ 
                         
                         ⁢ 
                         
                           
                             Z 
                             T 
                           
                           ⁡ 
                           
                             ( 
                             
                               
                                 
                                   
                                     E 
                                     T 
                                   
                                 
                                 
                                   
                                     E 
                                     0 
                                     T 
                                   
                                 
                               
                             
                             ) 
                           
                         
                       
                     
                     , 
                     
                       
                         y 
                         R 
                       
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 
                                   
                                     ( 
                                     
                                       
                                         A 
                                         1 
                                       
                                       ⁢ 
                                       
                                         R 
                                         1 
                                       
                                     
                                     ) 
                                   
                                   
                                     - 
                                     1 
                                   
                                 
                                 ⁢ 
                                 
                                   b 
                                   1 
                                 
                               
                             
                           
                           
                             
                               ⋮ 
                             
                           
                           
                             
                               
                                 
                                   
                                     ( 
                                     
                                       
                                         A 
                                         
                                           N 
                                           r 
                                         
                                       
                                       ⁢ 
                                       
                                         R 
                                         
                                           N 
                                           r 
                                         
                                       
                                     
                                     ) 
                                   
                                   
                                     - 
                                     1 
                                   
                                 
                                 ⁢ 
                                 
                                   b 
                                   
                                     N 
                                     r 
                                   
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       b 
                       ~ 
                     
                     = 
                     
                       
                         - 
                         
                           ( 
                           
                             
                               
                                 
                                   E 
                                   ⁢ 
                                   
                                     y 
                                     ^ 
                                   
                                 
                               
                             
                             
                               
                                 0 
                               
                             
                           
                           ) 
                         
                       
                       + 
                       
                         
                           ( 
                           
                             
                               
                                 E 
                               
                             
                             
                               
                                 
                                   E 
                                   0 
                                 
                               
                             
                           
                           ) 
                         
                         ⁢ 
                         
                           
                             ( 
                             
                               
                                 Ry 
                                 R 
                               
                               - 
                               
                                 Z 
                                 ⁢ 
                                 
                                   Q 
                                   ~ 
                                 
                                 ⁢ 
                                 
                                   
                                     Z 
                                     T 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     
                                       q 
                                       + 
                                       
                                         QRy 
                                         R 
                                       
                                     
                                     ) 
                                   
                                 
                               
                             
                             ) 
                           
                           . 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   23 
                   ) 
                 
               
             
           
         
       
     
     The full solution y cg ,λ cg  can be obtained as 
     
       
         
           
             
               
                 
                   
                     
                       y 
                       cg 
                     
                     = 
                     
                       
                         Ry 
                         R 
                       
                       - 
                       
                         Z 
                         ⁢ 
                         
                           Q 
                           ~ 
                         
                         ⁢ 
                         
                           
                             Z 
                             T 
                           
                           ⁡ 
                           
                             ( 
                             
                               
                                 
                                   
                                     E 
                                     T 
                                   
                                 
                                 
                                   
                                     E 
                                     0 
                                     T 
                                   
                                 
                               
                             
                             ) 
                           
                         
                       
                       - 
                       
                         Z 
                         ⁢ 
                         
                           Q 
                           ~ 
                         
                         ⁢ 
                         
                           
                             Z 
                             T 
                           
                           ⁡ 
                           
                             ( 
                             
                               q 
                               + 
                               
                                 QRy 
                                 R 
                               
                             
                             ) 
                           
                         
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       
                         λ 
                         1 
                         cg 
                       
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 
                                   E 
                                   1 
                                   T 
                                 
                                 ⁢ 
                                 
                                   
                                     λ 
                                     ^ 
                                   
                                   1 
                                 
                               
                             
                           
                           
                             
                               
                                 
                                   ∑ 
                                   2 
                                   
                                     N 
                                     r 
                                   
                                 
                                 ⁢ 
                                 
                                   
                                     λ 
                                     ~ 
                                   
                                   r 
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                     , 
                     
                       
                         λ 
                         r 
                         cg 
                       
                       = 
                       
                         
                           
                             ( 
                             
                               
                                 
                                   
                                     
                                       E 
                                       r 
                                       T 
                                     
                                     ⁢ 
                                     
                                       
                                         λ 
                                         ^ 
                                       
                                       r 
                                     
                                   
                                 
                               
                               
                                 
                                   
                                     
                                       λ 
                                       ~ 
                                     
                                     r 
                                   
                                 
                               
                             
                             ) 
                           
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             ∀ 
                             r 
                           
                         
                         = 
                         2 
                       
                     
                     , 
                     … 
                     ⁢ 
                     
                         
                     
                     , 
                     
                       
                         N 
                         r 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   24 
                   ) 
                 
               
             
           
         
       
     
     Check for Violation of Multiplier Sign 
     The solution (y cg ,λ cg ) is the optimal of Equation (6). Given a solution from the CG method, the solution is checked  720  to determine if all multipliers are of the correct sign by. The violated indices are found by
 
   I     r   viol   ={i∈ I     r |λ r,i &lt;0}, Ī   r   viol   ={i∈Ī   r |λ r,i &gt;0}.
 
     If the set  I   viol ∪Ī viol ≠Ø, then the active index sets  I   r ,Ī r  are updated  722  as,
 
   I     r   = I     r   \ I     r   viol   ,Ī   r   =Ī   r   \Ī   r   viol .
 
     If no such multiplier indices are found, then the obtained CG solution is checked  725  to see if it makes progress toward solving Equation (6). If sufficient progress is not made then, then the CG solution is discarded  730  and the procedure stops  750  and returns to the ADMM restarting from the previous ADMM iteration value. 
     Otherwise, check  735  of the termination condition is satisfied, and if yes, use  745  the CG solution to restart the ADMM and go to step  750 , and otherwise, select  740  the CG solution as a starting guess and continue the CG iterations. Namely, the iterates with superscripts  cg,1  and  cg,2  are used in  151  and  171 . 
     To do this, two iterations of the ADMM procedure are performed as follows,
 
 w   cg =             Y ( y   cg −λ cg /β)
 
( y   cg,1   ,w   cg,1 ,λ cg,1 )=         ( y   cg   ,w   cg ,λ cg /β)
 
( y   cg,2   ,w   cg,2 ,λ cg,2 )=         ( y   cg,1   ,w   cg,1 ,λ cg,1 )′
 
 v   cg,1   =y   cg,1 −λ cg   ,v   cg,2   =y   cg,2 −λ cg,1  

     Check for Progress Towards Solution of StQP 
     The condition checked for sufficient progress is,
 
∥ v   cg,2   −v   cg,1 ∥≦(1−η)∥ v   l   −v   l−1 ∥,
 
where η&lt;&lt;1 is a constant. If sufficient progress is made towards an optimal solution then, ADMM iterates are updated as,
 
( y   l   ,w   l ,λ l )=( y   cg,2   ,w   cg,2 ,λ cg,2 )
 
( y   l   ,w   l ,λ l )=( y   cg,2   ,w   cg,2 ,λ cg,2 )
 
 v   l   =v   cg,2   ,v   l−1   =v   cg,1 ,λ=βλ cg,2 .
 
     Namely, the iterates with superscripts  cg,1  and  cg,2  are used in  151  and  171 . If the termination conditions are satisfied then the ADDM-CG process  700  is terminated. If not, more CG iterations are performed using the computed CG iterates as initial solution  740 . The procedure is designed so that after the correct set of active indices are found, then the ADDM-CG process is used to compute the solution to Equation (6) and ADMM iterations are not used. 
     In the Appendix, pseudocode for the CG-based ADMM procedure is provided as Algorithm 1. Pseudocode for the CG process is provided in Algorithm 2, and pseudocode for finding indices of multipliers with the incorrect sign is provided as Algorithm 3. 
     Although the invention has been described by way of examples of preferred embodiments, it is to be understood that various other adaptations and modifications can be made within the spirit and scope of the invention. Therefore, it is the object of the appended claims to cover all such variations and modifications as come within the true spirit and scope of the invention. 
     
       
         
               
             
               
               
             
               
               
             
               
               
               
             
               
               
               
               
             
               
               
               
             
               
               
               
               
             
           
               
                 APPENDIX 
               
               
                   
               
               
                 Algorithm 1: (y l , w l , λ l  = ADMM-CG(y l , w l , λ l , v l , v l−1 ) 
               
               
                   
               
             
             
               
                   
               
             
          
           
               
                 Data:  
                 ADMM iterate (y l , w l , λ l ), v l  and v l-1   
               
               
                   
                 Active-set identification tolerance ε act ,  
               
               
                   
                 Convergence tolerance ε 
               
               
                   
                 Maximum number of CG iterations n cg, max , 
               
               
                   
                 Sufficient decrease constant η &lt;&lt; 1. 
               
               
                   
                 ADMM parameter β, 
               
             
          
           
               
                  1 
                 Set j = 0, flag = 0, resolvecg = 1 
               
               
                  2 
                 Define the active index sets as in (21) 
               
               
                  3 
                 Set λ = βλ 1   
               
               
                  4 
                 while flag == 0 do 
               
             
          
           
               
                  5 
                 | 
                 while resolvecg == 1 do 
               
             
          
           
               
                  6 
                 | 
                 | 
                 Call (y cg , λ cg ) = CG ( I , Ī, λ, n cg, max , ε) 
               
               
                  7 
                 | 
                 | 
                 Call ( I   viol , Ī viol ) = FindViolations(y cg , λ cg ) 
               
               
                  8 
                 | 
                 | 
                 if ( I   viol  = ∅ and Ī viol  = ∅) then 
               
               
                  9 
                 | 
                 | 
                 | Set resolvecg == 0 
               
               
                 10 
                 | 
                 | 
                 else 
               
               
                 11 
                 | 
                 | 
                 | Set  I   r  =  I   r \ I   r   viol , Ī r  = Ī r \Ī r   viol   
               
               
                 12 
                 | 
                 | 
                 end 
               
             
          
           
               
                 13 
                 | 
                 end 
               
               
                 14 
                 | 
                 Set w cg  = P Υ  (y cg  − λ cg /β) 
               
               
                 15 
                 | 
                 (y cg,1 , w cg,l , λ cg,l ) = ADMM(y cg , w cg , λ cg /β) 
               
               
                 16 
                 | 
                 (y cg,1 , w cg,2 ,λ cg,2 ) = ADMM(y cg,1  w cg,1 − λ cg,1 ) 
               
               
                 17 
                 | 
                 Set v cg,1 = y cg,1  − λ cg , v cg,2  − λ cg,1   
               
               
                 18 
                 | 
                 if ∥v cg,2  − v cg,1  ∥≦ (1-η)∥ v l -v l−1 ∥ then 
               
             
          
           
               
                 19 
                 | 
                 |  
                 Set (y l , w l , λ l  ) = (y cg,2 , w cg,2 , λ cg,2 ) 
               
               
                 20 
                 | 
                 | 
                 Set (y l , w l , λ l ) = (y cg,2  = w cg,2 , λ cg,2 ) 
               
               
                 21 
                 | 
                 | 
                 Set v l  = vc g,2 , v l−1  = v cg,1 , λ = βλ cg,2   
               
               
                 22 
                 | 
                 | 
                 if Termination conditions satisfied with (y cg,1 , w cg,1 , λ cg,1 ), 
               
               
                   
                 | 
                 | 
                 (y cg,2 , w cg,2 , λ cg,2 ) then 
               
               
                 23 
                 | 
                 | 
                 | flag = 1 
               
               
                 24 
                 | 
                 | 
                 end 
               
               
                 25 
                 | 
                 else 
                   
               
               
                 26 
                 | 
                 | 
                 flag = 1 
               
               
                 27 
                 | 
                 end 
                   
               
               
                 28 
                 end 
               
               
                   
               
             
          
         
       
     
     
       
         
               
             
               
               
               
               
             
               
               
               
               
               
             
               
               
               
               
             
           
               
                   
               
               
                 Algorithm 2: (y cg , λ cg ) = CG( I , Ī, λ, n cg,max , ε) 
               
               
                   
               
             
             
               
                   
               
             
          
           
               
                   
                   
                  1 
                 Form {tilde over (M)}, {tilde over (b)} as given in (23) 
               
               
                   
                   
               
               
                   
                   
                  2 
                 
                   
                     
                       
                         
                           Set 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             x 
                             0 
                           
                         
                         = 
                         
                           
                             ( 
                             
                               
                                 
                                   E 
                                 
                               
                               
                                 
                                   
                                     E 
                                     0 
                                   
                                 
                               
                             
                             ) 
                           
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           λ 
                         
                       
                     
                   
                 
               
               
                   
                   
               
               
                   
                   
                  3 
                 Set r 0  := {tilde over (b)} − {tilde over (M)}x 0 , p 0  := r 0 , k := 0 
               
               
                   
                   
                  4 
                 while ∥r k ∥ &gt; ε and k &lt; n cg,max  do 
               
               
                   
                   
               
             
          
           
               
                   
                   
                  5 
                 | 
                 
                   
                     
                       
                         
                           Set 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             α 
                             k 
                           
                         
                         ⁢ 
                         
                             
                         
                         := 
                         
                           
                             
                               r 
                               k 
                               T 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               r 
                               k 
                             
                           
                           
                             
                               p 
                               k 
                               T 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               Mp 
                               k 
                             
                           
                         
                       
                     
                   
                 
               
               
                   
                   
               
               
                   
                   
                  6 
                 | 
                 Set x k+1  := x k  + α k p k   
               
               
                   
                   
                  7 
                 | 
                 Set r k+1  := r k  − α k Mp k   
               
               
                   
                   
               
               
                   
                   
                  8 
                 | 
                 
                   
                     
                       
                         
                           Set 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             β 
                             k 
                           
                         
                         ⁢ 
                         
                             
                         
                         := 
                         
                           
                             
                               r 
                               
                                 k 
                                 + 
                                 1 
                               
                               T 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               r 
                               
                                 k 
                                 + 
                                 1 
                               
                             
                           
                           
                             
                               r 
                               k 
                               T 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               r 
                               k 
                             
                           
                         
                       
                     
                   
                 
               
               
                   
                   
               
               
                   
                   
                  9 
                 | 
                 Set p k+1  := r k+1  + β k p k   
               
               
                   
                   
                 10 
                 | 
                 Set k := k + 1 
               
             
          
           
               
                   
                   
                 11 
                 end 
               
               
                   
                   
                 12 
                 Compute (y cg , λ cg ) using (24). 
               
               
                   
                   
               
             
          
         
       
     
     
       
         
               
             
               
               
               
             
               
               
               
               
             
           
               
                   
               
               
                 Algorithm 3: ( I   viol , Ī viol ) = FindViolations(y, λ) 
               
               
                   
               
             
             
               
                   
               
             
          
           
               
                   
                 1 
                 for r = 1, . . . , N r  do 
               
             
          
           
               
                   
                 2 
                 | 
                 Set  I   r   viol  = {i ε  I   r  | λ r,i  &lt; 0} 
               
               
                   
                 3 
                 | 
                 Set Ī r   viol  = {i ε Ī r |  λ r,i  &gt; 0} 
               
               
                   
                 4 
                 end