Optimal parameter selection and acceleration in ADMM for multi-stage stochastic convex quadratic programs

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.

Claims

1. A method for controlling a machine by solving a stochastic quadratic program (StQP), wherein the StQP enumerates scenarios, wherein the StQP is multi-stage, stochastic and convex, wherein constraints of the StQP include a set of linear equalities for each scenario, a set of linear inequalities for each scenario, a set of non-anticipativity constraints for the scenarios, wherein the solving uses non-anticipativity constrained variables for each scenario and an Alternating Direction Method of Multipliers (ADMM) wherein the ADMM is applied to a model predictive controller (MPC) for a machine governed by a dynamical system, wherein variables of the StQP indicating desired behavior of the machine include a linear subspace constrained variable vector and a set constrained variable vector, wherein constraints include the set of linear inequalities in each scenario and the set of non-anticipativity constraints on the non-anticipativity constrained variables, comprising: solving, using a processor, the linear subspace constrained variable vector while keeping the set constrained variable vector fixed using an optimal step size for each scenario and a Lagrangian multiplier; solving, using the processor, the set constrained variable vector while keeping linear subspace constrained variable vector fixed using the optimal step size and the Lagrangian multiplier, wherein the non-anticipativity constrained variables satisfies the set of non-anticipativity constraints for the scenarios; updating the Lagrangian multiplier; repeating for a next iteration of the StQP until either a termination condition for a feasible solution is satisfied or a termination condition for an infeasible solution is satisfied to produce a StQP solution, wherein the StQP uses the ADMM with a conjugate gradient (CG) process (ADMM-CG), wherein the ADMM-CG is called after a predetermined number of iterations, further comprising: identifying a subset of linear inequalities that hold as equalities at a solution; performing the predetermined number of iterations of the CG process on a CG-linear system; performing 2 iterations of ADMM using the solution identified by CG process to determine if sufficient progress is made toward solving StQP; if sufficient progress is achieved, continuing with additional CG interations, and if sufficient progress is not achieved, then terminating the CG process; and controlling the machine using the StQP solution.

2. The method in claim 1 where the StQP is solved for each scenario, wherein the dynamical system depends on uncertain parameters.

3. The method of in claim 1 where the StQP uses the ADMM with a conjugate gradient (CG) process (ADMM-CG), wherein the ADMM-CG is called after a predetermined number of iterations, further comprising: identifying the subset of linear inequalities in the convex set that are expected hold as equalities at a solution; performing the predetermined number of iterations of the CG process on the CG-linear system; checking if all the multipliers for the identified set of linear inequalities are of the correct sign; if there exists multipliers that are not of the correct sign, then such linear inequalities are removed from the step and the CG iterations are repeated until all multipliers in such identified are of a correct sign or the set of identified linear inequalities is empty; performing 2 iterations of ADMM using the solution identified by CG process to determine if sufficient progress is made toward solving StQP; if sufficient progress is achieved, continuing with additional CG interations, and if sufficient progress is not achieved, then terminating the CG process.

4. The method in claim 1, wherein the optimal step size for each scenario is determined as a square root of a product of minimum and maximum eigenvalues of a Hessian matrix of the scenario pre- and post-multiplied by an orthonormal basis for a null space of a linear equality constraint matrix corresponding to the scenario.

5. The method in claim 1, wherein the optimal step size is determined as a square root of a product of minimum and maximum eigenvalues of a Hessian matrix of a QP problem pre- and post-multiplied by an orthonormal basis for a null space of a linear equality constraint matrix.

6. The method of claim 1, wherein the solving for the linear subspace constrained variable vector minimizes a Lagrangian augmented quadratic objective function depending on the Lagrange multiplier, and a convex set is subject to the set of linear equality constraints and the set of inequality constraints.

7. The method of claim 1, wherein the solving for the set constrained variable vector minimizes a Lagrangian augmented quadratic objective function depending on the Lagrange multiplier and the linear subspace constrained variable vector subject to the set of linear equality constraints and the set of inequality constraints.

8. The method of claim 1, wherein the updating uses a difference between the linear subspace constrained variable vector and the set constrained variable vector multiplied by an optimal step size parameter.

9. The method of claim 1, wherein the termination condition for the feasible solution satisfies a tolerance greater than zero of a maximum of a norm of a change in the set constrained variable vector from a current value to a value at a previous iteration, and a norm of a change in the Lagrange multiplier from a current value to a value at a previous iteration.

10. The method of claim 1, wherein the termination condition for the infeasible solution checks: satisfaction to a tolerance greater than zero of a maximum of a normed change in the set constrained variable vector from a current value to a value at a previous iteration; satisfaction to a tolerance greater than zero of the normed change in the linear subspace constrained variable vector from the current value to a value at the previous iteration; satisfaction to a tolerance greater than 0 of an angle between a Lagrange multiplier vector and a vector resulting from a difference of the linear subspace constrained variable vector and set the constrained variable vector at the current value; and 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.

11. The method of claim 1, wherein the linear subspace constrained variable vector includes states and controls of the dynamical system, and, optionally, outputs and constrained outputs along a prediction interval, and a set constrained auxiliary variable vector of a same size as the linear subspace constrained variable vector, wherein the linear subspace constraints are equations of prediction model dynamics along the prediction interval, and the set constraints are upper bounds and lower bounds on the states and the controls, and, optionally, upper bounds and lower bounds on the output and the constrained outputs along the prediction interval.

12. The method of claim 1, wherein a number of scenarios in the StQP is one and a number of non-anticipativity constraints is zero.

13. A method for controlling a machine by solving a stochastic quadratic program (StQP), wherein the StQP enumerates scenarios, wherein the StQP is multi-stage, stochastic and convex, wherein constraints of the StQP include a set of linear equalities for each scenario, a set of linear inequalities for each scenario, a set of non-anticipativity constraints for the scenarios, wherein the solving uses non-anticipativity constrained variables for each scenario and an Alternating Direction Method of Multipliers (ADMM) wherein the ADMM is applied to a model predictive controller (MPC) for a machine governed by a dynamical system, wherein variables of the StQP indicating desired behavior of the machine include a linear subspace constrained variable vector and a set constrained variable vector, wherein constraints include the set of linear inequalities in each scenario and the set of non-anticipativity constraints on the non-anticipativity constrained variables, comprising iterative steps: solving, using a processor, the linear subspace constrained variable vector while keeping the set constrained variable vector fixed using an optimal step size for each scenario and a Lagrangian multiplier; solving, using the processor, the set constrained variable vector while keeping linear subspace constrained variable vector fixed using the optimal step size and the Lagrangian multiplier, wherein the non-anticipativity constrained variables satisfies the set of non-anticipativity constraints for the scenarios; updating the Lagrangian multiplier; repeating for a next iteration of the StQP until either a termination condition for a feasible solution is satisfied or a termination condition for an infeasible solution is satisfied to produce a StQP solution; and controlling the machine using the StQP solution, where the StQP uses the ADMM with a conjugate gradient (CG) process (ADMM-CG), wherein the ADMM-CG is called after a predetermined number of iterations, further comprising: identifying the subset of linear inequalities in the convex set that are expected hold as equalities at a solution; performing the predetermined number of iterations of the CG process on the CG-linear system; checking if all the multipliers for the identified set of linear inequalities are of the correct sign; if there exists multipliers that are not of the correct sign, then such linear inequalities are removed from the step and the CG iterations are repeated until all multipliers in such identified are of a correct sign or the set of identified linear inequalities is empty; performing 2 iterations of ADMM using the solution identified by CG process to determine if sufficient progress is made toward solving StQP; if sufficient progress is achieved, continuing with additional CG interations, and if sufficient progress is not achieved, then terminating the CG process.

14. A method for controlling a machine by solving a stochastic quadratic program (StQP), wherein the StQP enumerates scenarios, wherein the StQP is multi-stage, stochastic and convex, wherein constraints of the StQP include a set of linear equalities for each scenario, a set of linear inequalities for each scenario, a set of non-anticipativity constraints for the scenarios, wherein the solving uses non-anticipativity constrained variables for each scenario and an Alternating Direction Method of Multipliers (ADMM) wherein the ADMM is applied to a model predictive controller (MPC) for a machine governed by a dynamical system, wherein variables of the StQP indicating desired behavior of the machine include a linear subspace constrained variable vector and a set constrained variable vector, wherein constraints include the set of linear inequalities in each scenario and the set of non-anticipativity constraints on the non-anticipativity constrained variables, comprising iterative steps: solving, using a processor, the linear subspace constrained variable vector while keeping the set constrained variable vector fixed using an optimal step size for each scenario and a Lagrangian multiplier; solving, using the processor, the set constrained variable vector while keeping linear subspace constrained variable vector fixed using the optimal step size and the Lagrangian multiplier, wherein the non-anticipativity constrained variables satisfies the set of non-anticipativity constraints for the scenarios; updating the Lagrangian multiplier; repeating for a next iteration of the StQP until either a termination condition for a feasible solution is satisfied or a termination condition for an infeasible solution is satisfied to produce a StQP solution, wherein the termination condition for the feasible solution satisfies a tolerance greater than zero of a maximum of a norm of a change in the set constrained variable vector from a current value to a value at a previous iteration, and a norm of a change in the Lagrange multiplier from a current value to a value at a previous iteration; and controlling the machine using the StQP solution.

15. A method for controlling a machine by solving a stochastic quadratic program (StQP), wherein the StQP enumerates scenarios, wherein the StQP is multi-stage, stochastic and convex, wherein constraints of the StQP include a set of linear equalities for each scenario, a set of linear inequalities for each scenario, a set of non-anticipativity constraints for the scenarios, wherein the solving uses non-anticipativity constrained variables for each scenario and an Alternating Direction Method of Multipliers (ADMM) wherein the ADMM is applied to a model predictive controller (MPC) for a machine governed by a dynamical system, wherein variables of the StQP indicating desired behavior of the machine include a linear subspace constrained variable vector and a set constrained variable vector, wherein constraints include the set of linear inequalities in each scenario and the set of non-anticipativity constraints on the non-anticipativity constrained variables, comprising iterative steps: solving, using a processor, the linear subspace constrained variable vector while keeping the set constrained variable vector fixed using an optimal step size for each scenario and a Lagrangian multiplier; solving, using the processor, the set constrained variable vector while keeping linear subspace constrained variable vector fixed using the optimal step size and the Lagrangian multiplier, wherein the non-anticipativity constrained variables satisfies the set of non-anticipativity constraints for the scenarios; updating the Lagrangian multiplier; repeating for a next iteration of the StQP until either a termination condition for a feasible solution is satisfied or a termination condition for an infeasible solution is satisfied to produce a StQP solution, wherein the termination condition for the infeasible solution checks: satisfaction to a tolerance greater than zero of a maximum of a normed change in the set constrained variable vector from a current value to a value at a previous iteration; satisfaction to a tolerance greater than zero of the normed change in the linear subspace constrained variable vector from the current value to a value at the previous iteration; satisfaction to a tolerance greater than 0 of an angle between a Lagrange multiplier vector and a vector resulting from a difference of the linear subspace constrained variable vector and set the constrained variable vector at the current value; and 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; and controlling the machine using the StQP solution.

16. A method for controlling a machine by solving a stochastic quadratic program (StQP), wherein the StQP enumerates scenarios, wherein the StQP is multi-stage, stochastic and convex, wherein constraints of the StQP include a set of linear equalities for each scenario, a set of linear inequalities for each scenario, a set of non-anticipativity constraints for the scenarios, wherein the solving uses non-anticipativity constrained variables for each scenario and an Alternating Direction Method of Multipliers (ADMM) wherein the ADMM is applied to a model predictive controller (MPC) for a machine governed by a dynamical system, wherein variables of the StQP indicating desired behavior of the machine include a linear subspace constrained variable vector and a set constrained variable vector, wherein constraints include the set of linear inequalities in each scenario and the set of non-anticipativity constraints on the non-anticipativity constrained variables, comprising iterative steps: solving, using a processor, the linear subspace constrained variable vector while keeping the set constrained variable vector fixed using an optimal step size for each scenario and a Lagrangian multiplier; solving, using the processor, the set constrained variable vector while keeping linear subspace constrained variable vector fixed using the optimal step size and the Lagrangian multiplier, wherein the non-anticipativity constrained variables satisfies the set of non-anticipativity constraints for the scenarios; updating the Lagrangian multiplier; repeating for a next iteration of the StQP until either a termination condition for a feasible solution is satisfied or a termination condition for an infeasible solution is satisfied to produce a StQP solution, wherein the linear subspace constrained variable vector includes states and controls of the dynamical system, and, optionally, outputs and constrained outputs along a prediction interval, and a set constrained auxiliary variable vector of a same size as the linear subspace constrained variable vector, wherein the linear subspace constraints are equations of prediction model dynamics along the prediction interval, and the set constraints are upper bounds and lower bounds on the states and the controls, and, optionally, upper bounds and lower bounds on the output and the constrained outputs along the prediction interval; and controlling the machine using the StQP solution.

17. A method for controlling a machine by solving a stochastic quadratic program (StQP), wherein the StQP enumerates scenarios, wherein the StQP is multi-stage, stochastic and convex, wherein constraints of the StQP include a set of linear equalities for each scenario, a set of linear inequalities for each scenario, a set of non-anticipativity constraints for the scenarios, wherein the solving uses non-anticipativity constrained variables for each scenario and an Alternating Direction Method of Multipliers (ADMM) wherein the ADMM is applied to a model predictive controller (MPC) for a machine governed by a dynamical system, wherein variables of the StQP indicating desired behavior of the machine include a linear subspace constrained variable vector and a set constrained variable vector, wherein constraints include the set of linear inequalities in each scenario and the set of non-anticipativity constraints on the non-anticipativity constrained variables, comprising iterative steps: solving, using a processor, the linear subspace constrained variable vector while keeping the set constrained variable vector fixed using an optimal step size for each scenario and a Lagrangian multiplier; solving, using the processor, the set constrained variable vector while keeping linear subspace constrained variable vector fixed using the optimal step size and the Lagrangian multiplier, wherein the non-anticipativity constrained variables satisfies the set of non-anticipativity constraints for the scenarios; updating the Lagrangian multiplier; repeating for a next iteration of the StQP until either a termination condition for a feasible solution is satisfied or a termination condition for an infeasible solution is satisfied to produce a StQP solution, wherein a number of scenarios in the StQP is one and a number of non-anticipativity constraints is zero; and controlling the machine using the StQP solution.

18. A method for controlling a machine by solving a stochastic quadratic program (StQP), wherein the StQP enumerates scenarios, wherein the StQP is multi-stage, stochastic and convex, wherein constraints of the StQP include a set of linear equalities for each scenario, a set of linear inequalities for each scenario, a set of non-anticipativity constraints for the scenarios, wherein the solving uses non-anticipativity constrained variables for each scenario and an Alternating Direction Method of Multipliers (ADMM) wherein the ADMM is applied to a model predictive controller (MPC) for the machine governed by a dynamical system, comprising: solving, using a processor, the StQP iteratively until either a termination condition for a feasible solution is satisfied or a termination condition for an infeasible solution is satisfied to produce a StQP solution including one or combination of values of the linear subspace constrained variable vector and the set of constrained variable vector, wherein variables of the StQP indicate desired behavior of the machine and include a linear subspace constrained variable vector subject to the constraints including the set of linear equalities for each scenario and a set constrained variable vector subject the constraints including the set of linear inequalities in each scenario for the variables excluding the non-anticipativity constrained variables and the set of non-anticipativity constraints on the non-anticipativity constrained variables, wherein each iteration comprises: solving the linear subspace constrained variable vector while keeping the set constrained variable vector fixed using an optimal step size for each scenario and a Lagrangian multiplier; solving the set constrained variable vector while keeping linear subspace constrained variable vector fixed using the optimal step size and the Lagrangian multiplier, wherein the non-anticipativity constrained variables satisfies the set of non-anticipativity constraints for the scenarios; and updating the Lagrangian multiplier; and controlling the machine using the StQP solution.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) 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;

(2) 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;

(3) FIG. 2B is a schematic of a convex region Y defined as hyperplanes that can be solved by embodiments of the invention;

(4) FIG. 3 is a block diagram of a controller and machine using the ADMM according to embodiments of the invention;

(5) FIG. 4 is a block diagram of the controller using the ADMM according to embodiments of the invention;

(6) FIGS. 5A and 5B are schematics of scenario generation for stochastic MPC control according to embodiments of the invention;

(7) FIG. 6 is a flow diagram of stochastic MPC based on scenario enumeration according to embodiments of the invention; and

(8) 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

(9) 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).

(10) 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.

(11) 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)

(12) 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.sub.m(w(t))x(t)+B.sub.m(w(t))u(t)+G.sub.m(w(t))
p(w(t))=ƒ.sub.p(w(t−1), . . . ,w(t−t.sub.w)),  (1)
where x∈.sup.n.sup.x is the difference Equation state which contains the machine state and information of the desired behavior of the machine, u∈.sup.n.sup.u is the input, w∈.sup.n.sup.w is a stochastic vector that represents the uncertainty, possibly with a finite discrete set of values allowed, p is a probability distribution function, ƒ.sub.p is a function that describes p, t.sub.w is a length of time steps t that defines the probability distribution, and A.sub.m, B.sub.m, G.sub.m are parameters of the model.

(13) The states and inputs of the system whose behavior is described in Equation (1) can be subject to constraints
x(t+k)∈X.sub.k|t,u(t+k)∈U.sub.k|t,  (2)
where X.sub.k|t, U.sub.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

(14) $\begin{array}{cc}\underset{u_t}{\min }E\left[x_{N|t}^T{\mathrm{Px}}_{N|t}+\underset{k=0}{\overset{N}{.Math.}}{x}_{k|t}^T{\mathrm{Qx}}_{k|t}+u_{k|t}^T{\mathrm{Ru}}_{k|t}\right]s.t.x_{k+1|t}=A_m\left(w\left(t+k\right)\right)x_{k|t}+B_m\left(w\left(t+k\right)\right)u_{k|t}+G_m\left(w\left(t+k\right)\right)p\left(w\left(t+k\right)\right)=f_p\left(w\left(t+k-1\right),•,w\left(t+k-t_w\right)\right)x_{k|t}\in X_{k|t},u_{k|t}\in U_{k|t}{x}_{0|t}=x\left(t\right),& \left(3\right)\end{array}$
where a.sub.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.sub.t=(u.sub.0|t, . . . , u.sub.N-1|t) is the control policy.

(15) 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

(16) 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

(17) The scenarios are sequences of possible realizations w, i.e.,
s.sub.k|t(j)=[w.sub.0|t.sup.j . . . w.sub.N-1|t.sup.j],
and to each scenario, a probability π(s.sub.k|t(j)) can be associated via p and ƒ.sub.p. The scenarios can be generated as described for FIGS. 5A and 5B. From a current estimate w(t) of the uncertainty 501 r.sub.1 possible values for the future uncertainty are selected w.sub.0|t.sup.j, j=1, . . . r.sub.1 for instance based on their likelihood as described by the function ƒ.sub.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

(18) 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.

(19) In scenario-enumeration MPC for N.sub.r scenarios, which is a subset of the possible scenarios on the MPC prediction horizon, the scenario-enumeration stochastic MPC problem is formulated as

(20) $\begin{array}{cc}\underset{{\left\{\left(x_t^r,u_t^r\right)\right\}}_{r=1}^{N_r}}{\min }\underset{r=0}{\overset{N_r}{.Math.}}\pi \left(s_{k|t}\left(r\right)\right)\left[{\left(x_{N|t}^r\right)}^T{\mathrm{Px}}_{N|t}^r+\underset{k=0}{\overset{N-1}{.Math.}}{\left(x_{k|t}^r\right)}^T{\mathrm{Qx}}_{k|t}^r+{\left(u_{k|t}^r\right)}^T{\mathrm{Ru}}_{k|t}^r\right]s.t.x_{k+1|t}=A_m\left(w_{k|t}^r\right)x_{k|t}^r+B_m\left(w_{k|t}^r\right)u_{k|t}^r+G_m\left(w_{k|t}^r\right)\forall r=1,\mathrm{.Math.},N_r,k=0,\mathrm{.Math.},N-1x_{k+1|t}^r\in X_{k+1|t},u_{k|t}^r\in U_{k|t}\forall r=1,\mathrm{.Math.},N_r,k=0,\mathrm{.Math.},N-1x_{0|t}^r=x\left(t\right)\forall r=1,\mathrm{.Math.},N_r{w}_{0|t}^r=w\left(t\right)\forall r=1,\mathrm{.Math.},N_r{u}_{k|t}^{r^{\min }\left(k\right)}=u_{k|t}^r,\forall r\ne r^{\min }\left(k\right),r\in R_{k|t}\mathrm{and}{r}^{\min }\left(k\right)=\min \left(R_{k|t}\right).& \left(4\right)\end{array}$
where x.sub.t.sup.r=(x.sub.1|t.sup.r . . . x.sub.N|t.sup.r), u.sub.t.sup.r=(u.sub.0|t.sup.r . . . u.sub.N-1|t.sup.r), and R.sub.k|t is the set of all scenarios, which are equal from the beginning until the step k−1, that is,
r,r′∈R.sub.k|ts.sub.k′|t(r)=s.sub.k′|t(r′)∀k′=0, . . . ,k−1.

(21) Due to the scenario enumeration, all vectors in Equation (4) are now deterministic, so that Equation (4) can be formulated as

(22) $\begin{array}{cc}\underset{{\left\{y_r\right\}}_{r=1}^{N_r}}{\min }\underset{r=1}{\overset{N_r}{.Math.}}\left(\frac{1}{2}{y}_r^T{Q}_r{y}_r+q_r^T{y}_r\right)s.t.A_r{y}_r=b_r\forall r=1,\mathrm{.Math.},N_r{y}_r\in Y_r\forall r=1,\mathrm{.Math.},N_r{\stackrel{\_}{y}}_k^{r^{\min }\left(k\right)}={\stackrel{\_}{y}}_k^r\forall r\ne r^{\min }\left(k\right),r\in R_{k|t}\mathrm{and}{r}^{\min }\left(k\right)=\min \left(R_{k|t}\right),& \left(5\right)\end{array}$
where for each scenario r=1, . . . , N.sub.r, k.sup.max(r) is the largest index for which r∈R.sub.k.sub.max.sub.(r)|t. The variables y.sub.k.sup.r are called the non-anticipativity constrained variables. Note that, k.sup.max(r)≧0, which means that at a minimum the control at the current time t must be equated across all scenarios and

(23) $y_r=\left(y^r,{\stackrel{\_}{y}}_0^r,\mathrm{.Math.},{\stackrel{\_}{y}}_k^r{\max }_{\left(r\right)}\right),y^r=\left(x_t^r=u_t^r\right),q_r=0$${\Pi }_x^r=\left(\begin{array}{ccc}\pi \left(s_{1|t}\left(r\right)\right)& & 0\\ & \ddots & \\ 0& \mathrm{.Math.}& \pi \left(s_{N-1|t}\left(r\right)\right)\end{array}\right)$${\Pi }_u^r=\left(\begin{array}{ccc}\pi \left(s_{0|t}\left(r\right)\right)& & 0\\ & \ddots & \\ 0& \mathrm{.Math.}& \pi \left(s_{N-1|t}\left(r\right)\right)\end{array}\right)$$Q_r=\left(\begin{array}{cccc}{\Pi }_x^r.Math.Q& 0& 0& 0\\ 0& \pi \left(s_{k|t}\left(r\right)\right)P& 0& 0\\ 0& 0& {\Pi }_u^r.Math.R& 0\\ 0& 0& 0& 0\end{array}\right)$$b=\left(\begin{array}{cc}{A}_m\left(w_{0|t}^r\right)x\left(t\right)+G_m\left(w_{0|t}^r\right)& \\ G_m\left(w_{1|t}^r\right)& \\ \mathrm{.Math.}& \\ G_m\left(w_{N-1|t}^r\right)& \\ \begin{array}{c}0\\ \mathrm{.Math.}\\ 0\end{array}\}& n_u\left(k^{\max }\left(r\right)+1\right)\mathrm{times}\end{array}\right)$$Y_r=\underset{\underset{=:Y_r}{︸}}{X_{1|t}\times \mathrm{.Math.}{X}_{N|t}\times U_{0|t}\times \mathrm{.Math.}{U}_{N-1|t}}\times R^{n_u\left(k^{\max }\left(r\right)+1\right)}.$

(24) The matrix A.sub.r is defined so that,

(25) $A_r=\left(\begin{array}{ccc}{\hat{A}}_r& {\hat{B}}_r& 0\\ 0& \left(\begin{array}{cc}-I_{n_u\left(k^{\max }\left(r\right)+1\right)}& 0\end{array}\right)& I_{n_u\left(k^{\max }\left(r\right)+1\right)}\end{array}\right)$$\mathrm{where},{\hat{A}}_r=\left(\begin{array}{ccccc}{I}_{n_x}& 0& \mathrm{.Math.}& 0& 0\\ -A_m\left(w_{1|t}^r\right)& I_{n_x}& \mathrm{.Math.}& 0& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& \ddots & \mathrm{.Math.}& \mathrm{.Math.}\\ 0& 0& \mathrm{.Math.}& I_{n_x}& 0\\ 0& 0& \mathrm{.Math.}& -A_m\left(w_{N-1|t}^r\right)& I_{n_x}\end{array}\right)$$\mathrm{and},{\hat{B}}_r=\left(\begin{array}{ccc}-B_m\left(w_{0|t}^r\right)& \mathrm{.Math.}& 0\\ & \ddots & \\ 0& \mathrm{.Math.}& -B_m\left(w_{N-1|t}^r\right)\end{array}\right),I_{n_u\left(k^{\max }\left(r\right)+1\right)}\mathrm{is}\mathrm{an}\mathrm{identity}\mathrm{matrix}\mathrm{of}\mathrm{size}{n}_u\left(k^{\max }\left(r\right)+1\right)\times n_u\left(k^{\max }\left(r\right)+1\right).$

Scenario Decomposition-Based ADMM

(26) 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.

(27) In one embodiment of the invention, the quadratic program (5) resulting from scenario-enumeration stochastic MPC problem can be reformulated as,

(28) $\begin{array}{cc}\underset{{\left\{y_r\right\}}_{r=1}^{N_r}}{\min }\underset{r=1}{\overset{N_r}{.Math.}}\left(\frac{1}{2}{y}_r^T{Q}_r{y}_r+q_r^T{y}_r\right)s.t.A_r{y}_r=b_r\forall r=1,\mathrm{.Math.},N_r,y\in Y_{\mathrm{sd}}& \left(6\right)\end{array}$
where a convex set Y.sub.sd is defined as

(29) $\begin{array}{cc}{Y}_{\mathrm{sd}}=\left\{\left(y_1,\mathrm{.Math.},y_{N_r}\right)|\begin{array}{c}{y}_r\in Y_r\forall r=1,\mathrm{.Math.},N_r\\ \begin{array}{c}{\stackrel{\_}{y}}_k^{r^{\min }}\left(k\right)={\stackrel{\_}{y}}^{r^{\min }\left(k\right)}\forall r\ne r_k,r\in R_{k|t}\\ \mathrm{and}{r}^{\min }\left(k\right)=\min \left(R_{k|t}\right)\end{array}\end{array}\right\}.& \left(7\right)\end{array}$

Convex Sets

(30) 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.sub.sd.

(31) Consider reformulating Equation (6) as,

(32) $\begin{array}{cc}\underset{y,w}{\min }\underset{r=1}{\overset{N_r}{.Math.}}\left(\frac{1}{2}{y}_r^T{Q}_r{y}_r+q_r^T{y}_r\right)s.t.A_r{y}_r=b_r\forall r=1,\mathrm{.Math.},N_r,y=w,w\in Y_{\mathrm{sd}}& \left(8\right)\end{array}$
where y=(y.sub.1, . . . , y.sub.N.sub.r)∈.sup.n with n=N.sub.r(n.sub.x+n.sub.u+n.sub.u(k.sup.max(r)+1)). The variables w∈.sup.n are required to be in the convex set Y.sub.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.

(33) The ADMM procedure dualizes the constraints y.sub.r=w.sub.r into the objective function using scaled multipliers β.sub.rλ.sub.r respectively where β.sub.r>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

(34) $\begin{array}{cc}\underset{y,w}{\min }\underset{r=1}{\overset{N_r}{.Math.}}{L}_r\left(y_r,w_r,{\lambda }_r\right)s.t.A_r{y}_r=b_r\forall r=1,\mathrm{.Math.},N_r,w\in Y_{\mathrm{sd}}& \left(9\right)\end{array}$
where

(35) 0$\begin{array}{cc}{L}_r\left(y_r,w_r,{\lambda }_r\right):=\frac{1}{2}{y}_r^T{Q}_r{y}_r+q_r^T{y}_r+\frac{{\beta }_r}{2}{.Math.y_r-w_r-{\lambda }_r.Math.}^2.& \left(10\right)\end{array}$

(36) 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:

(37) $\begin{array}{cc}{y}^{l+1}=\arg \underset{y}{\min }\underset{r=1}{\overset{N_r}{.Math.}}{L}_r\left(y_r,w_r^l,{\lambda }_r^l\right)s.t.A_r{y}_r=b_r\forall r=1,\mathrm{.Math.},N_r& \left(11\right)\\ w^{l+1}=\arg \underset{w}{\min }\underset{r=1}{\overset{N_r}{.Math.}}{L}_r\left(y_r^{l+1},w,{\tilde{\lambda }}_r^l\right)s.t.w\in Y_{\mathrm{sd}}& \left(12\right)\\ {\lambda }^{l+1}={\lambda }^l+w^{l+1}-y^{l+1}.& \left(13\right)\end{array}$

(38) 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.

(39) Observe that the objective functions and constraints in Equation (11) are decoupled by scenarios. Hence, the update 111 can be rewritten as,
y.sub.r.sup.l+1=M.sub.r(w.sub.r.sup.l+λ.sub.r.sup.l−q.sub.r/β.sub.r)+N.sub.rb.sub.r
where,
M.sub.r:=Z.sub.r(Z.sub.r.sup.T(Q.sub.r/β.sub.r+I.sub.n)Z.sub.r).sup.−1Z.sub.r.sup.T,
N.sub.r:=(I.sub.n−M.sub.rQ.sub.r/β.sub.r)R.sub.r(A.sub.rR.sub.r).sup.−1,  (14)
with R.sub.r,Z.sub.r denote an orghonormal basis for the range space of A.sub.r.sup.T and null space of A.sub.r, respectively. Thus, the update for y decouples by scenario and scales linearly with the number of scenarios.

(40) In (12), the objective function is component wise separable in the w. Using w.sub.r=(w.sub.r,w.sub.r) the constraint set w∈Y.sub.sd can be posed as,
w.sup.r∈Y.sub.r∀r=1, . . . ,n.sub.r and
w.sub.k.sup.r.sup.min.sup.(k)=w.sub.k.sup.r∀r≠r.sup.min(k),r∈R.sub.k|t and r.sup.min(k)=min(R.sub.k|t)  (15)

(41) Note that w.sub.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

(42) $\begin{array}{cc}{w}_r^{l+1}=P_Y\left({\left(y^r\right)}^{l+1}-{\left({\lambda }^r\right)}^l\right){\left({\stackrel{\_}{w}}_k^r\right)}^{l+1}=\frac{1}{.Math.R_{k|t}.Math.}\underset{r^{\prime }\in R_{k|t}}{.Math.}\left({\left({\stackrel{\_}{y}}_k^{r^{\prime }}\right)}^{l+1}-{\left({\stackrel{\_}{\lambda }}^r\right)}^l\right)\forall k=,0,\mathrm{.Math.},N-1,r\in R_{k|t}& \left(16\right)\end{array}$
where we have used λ.sub.r=(λ.sup.r,λ.sup.r) and P.sub.Y(x) denotes the projection of vector x onto the set Y. Even though the set Y.sub.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

(43) A check is made 133 to call the CG process every n.sup.admm iterations, where n.sup.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

(44) Step 150 checks if an optimality termination holds. Given a predetermined tolerance ∈, a termination condition 151 for an optimal feasible solution 160 is
r.sup.opt=max(∥w.sup.l+1−w.sup.l∥,∥λ.sup.l+1−λ.sup.l∥)<∈.  (17)

(45) 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.

(46) Otherwise, for ∈>0, a termination condition 171 for certifying a solution to the problem is infeasibility is checked in step 170

(47) $\begin{array}{cc}{r}^{\inf }=\max \left(.Math.y^{l+1}-y^l.Math.,.Math.w^{l+1}-w^l.Math.,1-\frac{{\left(w^{l+1}-y^{l+1}\right)}^T{\lambda }^{l+1}}{.Math.w^{l+1}-y^{l+1}.Math..Math.{\lambda }^{l+1}.Math.}\right)<\mathrm{.Math.}\left(y^{l+1}-w^{l+1}\right)\circ {\lambda }^{l+1}\ge 0& \left(18\right)\end{array}$
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.

(48) The termination condition 171 in Equation (18) is checked for the satisfaction of four conditions.

(49) 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.

(50) 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.

(51) 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.

(52) 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.

(53) 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

(54) The choice of the step size β.sub.r, which ensures that a least number of iterations are required for termination of the ADMM method is
β.sub.r.sup.opt=√{square root over (λ.sub.min(Z.sub.r.sup.TQ.sub.rZ.sub.r)λ.sub.max(Z.sub.r.sup.TQ.sub.rZ.sub.r))},  (19)
where λ.sub.min(•),λ.sub.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.

(55) In another embodiment of the invention, a single value of β is chosen for all scenarios β.sub.r=β. In this case, the optimal value of the parameter is prescribed as,
β.sup.opt=√{square root over (λ.sub.min(Z.sup.TQZ)λ.sub.max(Z.sup.TQZ))},  (20)
where, Z is the orthonormal basis for the vectors satisfying the constraints,
A.sub.ry.sub.r=0∀r=1, . . . ,N.sub.r
y.sub.k.sup.r.sup.min.sup.(k)=y.sub.k.sup.r∀r≠r.sub.k,r∈R.sub.k|t and r.sup.min(k)=min(R.sub.k|t)
and,

(56) $Q=\left(\begin{array}{ccc}{Q}_1& & 0\\ & \ddots & \\ 0& \mathrm{.Math.}& Q_{N_r}\end{array}\right).$

Scenario Decomposition and Conjugate Gradient (CG)-Based ADMM

(57) 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.sub.r=[y.sup.r,y.sup.r], where y.sup.r is a vector of lower bounds and y.sup.r is vector of upper bounds for y.sup.r, and the scenarios are such that R.sub.k|t=Ø∀k>1. In other words, the MPC problem is a two-stage stochastic program.

(58) Hence,
y.sub.r=(y.sup.r,y.sub.0.sup.r),w.sub.r=(w.sup.r,w.sub.0.sup.r),λ.sub.r=(λ.sup.r,λ.sub.0.sup.r).Math.r=1, . . . ,N.sub.r,
where
(y.sup.+,w.sup.+,λ.sup.+)=ADMM(y,w,λ)
denotes one iteration of ADMM, that is application of Equations (11-13), where (y.sup.k, w.sup.k, λ.sup.k)=(y, w, λ), and the parameter β.sub.r=β. For convenience, let v.sup.l=y.sup.l−λ.sup.l−1.

(59) After every n.sup.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

(60) Given an ADMM iterate (y.sup.l, w.sup.l, λ.sup.l) 705, we define the index sets 710 as

(61) $\begin{array}{cc}{\underset{\_}{I}}_r=\left\{i|{\left(w_i^r\right)}^l={\underset{\_}{y}}_i^r\mathrm{or}{\left({\lambda }_i^r\right)}^l\ge {\epsilon }^{\mathrm{act}}\right\}{\stackrel{\_}{I}}_r=\left\{i|{\left(w_i^r\right)}^l={\stackrel{\_}{y}}_i^r\mathrm{or}{\left({\lambda }_i^r\right)}^l\le -{\epsilon }^{\mathrm{act}}\right\}{\hat{y}}_{r,i}=\{\begin{array}{c}{\underset{\_}{y}}_i^r\mathrm{if}i\in {\underset{\_}{I}}_r\\ {\stackrel{\_}{y}}_i^r\mathrm{if}i\in {\stackrel{\_}{I}}_r\\ 0\mathrm{otherwise}\end{array},& \left(21\right)\end{array}$
where ∈.sup.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.sub.r is a matrix that extracts the components of w.sup.r that are in I.sub.r∪Ī.sub.r, and E.sub.r0 is a matrix that extracts the components of w.sub.r that correspond to y.sub.0.sup.r. The above active index sets in I.sub.r,Ī.sub.r are used as estimates of the set of indices of y.sub.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:

(62) $\begin{array}{cc}\tilde{M}\left(\begin{array}{c}\hat{\lambda }\\ \stackrel{\_}{\lambda }\end{array}\right)=\tilde{b}& \left(22\right)\end{array}$
where

(63) $\begin{array}{cc}\hat{\lambda }=\left(\begin{array}{c}{\hat{\lambda }}^1\\ \mathrm{.Math.}\\ {\hat{\lambda }}^{N_r}\end{array}\right),\hat{\lambda }=\left(\begin{array}{c}{\tilde{\lambda }}^2\\ \mathrm{.Math.}\\ {\tilde{\lambda }}^{N_r}\end{array}\right),q=\left(\begin{array}{c}{q}_1\\ \mathrm{.Math.}\\ q_{N_r}\end{array}\right),b=\left(\begin{array}{c}{b}_1\\ \mathrm{.Math.}\\ b_{N_r}\end{array}\right),Q=\left(\begin{array}{ccc}{Q}_1& & \\ & \ddots & \\ & & Q_{N_r}\end{array}\right),A=\left(\begin{array}{ccc}{A}_1& & \\ & \ddots & \\ & & A_{N_r}\end{array}\right),E=\left(\begin{array}{ccc}{E}_1& & \\ & \ddots & \\ & & E_{N_r}\end{array}\right),E_0=\left(\begin{array}{cccc}{E}_{10}& -E_{20}& 0& \\ E_{10}& 0& -E_{30}& \\ \mathrm{.Math.}& & & \ddots \\ E_{10}& & 0& -E_{n_{r}0}\end{array}\right)Z=\left(\begin{array}{ccc}{Z}_1& & \\ & \ddots & \\ & & Z_{N_r}\end{array}\right),R=\left(\begin{array}{ccc}{R}_1& & \\ & \ddots & \\ & & R_{N_r}\end{array}\right),\tilde{Q}=\left(\begin{array}{ccc}{\left(Z_1^T{Q}_r{Z}_1\right)}^{-1}& & \\ & \ddots & \\ & & {\left(Z_{N_r}^T{Q}_{N_r}{Z}_{N_r}\right)}^{-1}\end{array}\right)\tilde{M}=\left(\begin{array}{c}E\\ E_0\end{array}\right)Z\tilde{Q}{Z}^T\left(\begin{array}{cc}{E}^T& E_0^T\end{array}\right),y_R=\left(\begin{array}{c}{\left(A_1{R}_1\right)}^{-1}{b}_1\\ \mathrm{.Math.}\\ {\left(A_{N_r}{R}_{N_r}\right)}^{-1}{b}_{N_r}\end{array}\right)\tilde{b}=-\left(\begin{array}{c}E\hat{y}\\ 0\end{array}\right)+\left(\begin{array}{c}E\\ E_0\end{array}\right)\left({\mathrm{Ry}}_R-Z\tilde{Q}{Z}^T\left(q+{\mathrm{QRy}}_R\right)\right).& \left(23\right)\end{array}$

(64) The full solution y.sup.cg,λ.sup.cg can be obtained as

(65) $\begin{array}{cc}{y}^{\mathrm{cg}}={\mathrm{Ry}}_R-Z\tilde{Q}{Z}^T\left(\begin{array}{cc}{E}^T& E_0^T\end{array}\right)-Z\tilde{Q}{Z}^T\left(q+{\mathrm{QRy}}_R\right){\lambda }_1^{\mathrm{cg}}=\left(\begin{array}{c}{E}_1^T{\hat{\lambda }}^1\\ \underset{2}{\overset{N_r}{.Math.}}{\tilde{\lambda }}^r\end{array}\right),{\lambda }_r^{\mathrm{cg}}=\left(\begin{array}{c}{E}_r^T{\hat{\lambda }}^r\\ {\tilde{\lambda }}^r\end{array}\right)\forall r=2,\mathrm{.Math.},N_r.& \left(24\right)\end{array}$

Check for Violation of Multiplier Sign

(66) The solution (y.sup.cg,λ.sup.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.sub.r.sup.viol={i∈I.sub.r|λ.sub.r,i<0},Ī.sub.r.sup.viol={i∈Ī.sub.r|λ.sub.r,i>0}.

(67) If the set I.sup.viol∪Ī.sup.viol≠Ø, then the active index sets I.sub.r,Ī.sub.r are updated 722 as,
I.sub.r=I.sub.r\I.sub.r.sup.viol,Ī.sub.r=Ī.sub.r\Ī.sub.r.sup.viol.

(68) 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.

(69) 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 .sup.cg,1 and .sup.cg,2 are used in 151 and 171.

(70) To do this, two iterations of the ADMM procedure are performed as follows,
w.sup.cg=.sub.Y(y.sup.cg−λ.sup.cg/β)
(y.sup.cg,1,w.sup.cg,1,λ.sup.cg,1)=(y.sup.cg,w.sup.cg,λ.sup.cg/β)
(y.sup.cg,2,w.sup.cg,2,λ.sup.cg,2)=(y.sup.cg,1,w.sup.cg,1,λ.sup.cg,1)′
v.sup.cg,1=y.sup.cg,1−λ.sup.cg,v.sup.cg,2=y.sup.cg,2−λ.sup.cg,1

Check for Progress Towards Solution of StQP

(71) The condition checked for sufficient progress is,
∥v.sup.cg,2−v.sup.cg,1∥≦(1−η)∥v.sup.l−v.sup.l−1∥,
where η<<1 is a constant. If sufficient progress is made towards an optimal solution then, ADMM iterates are updated as,
(y.sup.l,w.sup.l,λ.sup.l)=(y.sup.cg,2,w.sup.cg,2,λ.sup.cg,2)
(y.sup.l,w.sup.l,λ.sup.l)=(y.sup.cg,2,w.sup.cg,2,λ.sup.cg,2)
v.sup.l=v.sup.cg,2,v.sup.l−1=v.sup.cg,1,λ=βλ.sup.cg,2.

(72) Namely, the iterates with superscripts .sup.cg,1 and .sup.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.

(73) 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.

(74) 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.

(75) TABLE-US-00001 APPENDIX Algorithm 1: (y.sup.l, w.sup.l, λ.sup.l = ADMM-CG(y.sup.l, w.sup.l, λ.sup.l, v.sup.l, v.sup.l−1) Data: ADMM iterate (y.sup.l, w.sup.l, λ.sup.l), v.sup.l and v.sup.l-1 Active-set identification tolerance ε.sup.act, Convergence tolerance ε Maximum number of CG iterations n.sup.cg, max, Sufficient decrease constant η << 1. ADMM parameter β,  1 Set j = 0, flag = 0, resolvecg = 1  2 Define the active index sets as in (21)  3 Set λ = βλ.sup.1  4 while flag == 0 do  5 | while resolvecg == 1 do  6 | | Call (y.sup.cg, λ.sup.cg) = CG (I, Ī, λ, n.sup.cg, max, ε)  7 | | Call (I.sup.viol, Ī.sup.viol) = FindViolations(y.sup.cg, λ.sup.cg)  8 | | if (I.sup.viol = ∅ and Ī.sup.viol = ∅) then  9 | | | Set resolvecg == 0 10 | | else 11 | | | Set I.sub.r = I.sub.r\I.sub.r.sup.viol, Ī.sub.r = Ī.sub.r\Ī.sub.r.sup.viol 12 | | end 13 | end 14 | Set w.sup.cg = P.sub.Υ (y.sup.cg − λ.sup.cg/β) 15 | (y.sup.cg,1, w.sup.cg,l, λ.sup.cg,l) = ADMM(y.sup.cg, w.sup.cg, λ.sup.cg/β) 16 | (y.sup.cg,1, w.sup.cg,2,λ.sup.cg,2) = ADMM(y.sup.cg,1 w.sup.cg,1− λ.sup.cg,1) 17 | Set v.sup.cg,1= y.sup.cg,1 − λ.sup.cg, v.sup.cg,2 − λ.sup.cg,1 18 | if ∥v.sup.cg,2 − v.sup.cg,1 ∥≦ (1-η)∥ v.sup.l-v.sup.l−1∥ then 19 | | Set (y.sup.l, w.sup.l, λ.sup.l ) = (y.sup.cg,2, w.sup.cg,2, λ.sup.cg,2) 20 | | Set (y.sup.l, w.sup.l, λ.sup.l) = (y.sup.cg,2 = w.sup.cg,2, λ.sup.cg,2) 21 | | Set v.sup.l = vc.sup.g,2, v.sup.l−1 = v.sup.cg,1, λ = βλ.sup.cg,2 22 | | if Termination conditions satisfied with (y.sup.cg,1, w.sup.cg,1, λ.sup.cg,1), | | (y.sup.cg,2, w.sup.cg,2, λ.sup.cg,2) then 23 | | | flag = 1 24 | | end 25 | else 26 | | flag = 1 27 | end 28 end

(76) TABLE-US-00002 Algorithm 2: (y.sup.cg, λ.sup.cg) = CG(I, Ī, λ, n.sup.cg,max, ε)    1 Form {tilde over (M)}, {tilde over (b)} as given in (23)  2 $\mathrm{Set}{x}_0=\left(\begin{array}{c}E\\ E_0\end{array}\right)\lambda$  3 Set r.sub.0 := {tilde over (b)} − {tilde over (M)}x.sub.0, p.sub.0 := r.sub.0, k := 0  4 while ∥r.sub.k∥ > ε and k < n.sup.cg,max do  5 | 0$\mathrm{Set}{\alpha }_k:=\frac{r_k^T{r}_k}{p_k^T{\mathrm{Mp}}_k}$  6 | Set x.sub.k+1 := x.sub.k + α.sub.kp.sub.k  7 | Set r.sub.k+1 := r.sub.k − α.sub.kMp.sub.k  8 | $\mathrm{Set}{\beta }_k:=\frac{r_{k+1}^T{r}_{k+1}}{r_k^T{r}_k}$  9 | Set p.sub.k+1 := r.sub.k+1 + β.sub.kp.sub.k 10 | Set k := k + 1 11 end 12 Compute (y.sup.cg, λ.sup.cg) using (24).

(77) TABLE-US-00003 Algorithm 3: (I.sup.viol, Ī.sup.viol) = FindViolations(y, λ) 1 for r = 1, . . . , N.sub.r do 2 | Set I.sub.r.sup.viol = {i ε I.sub.r | λ.sub.r,i < 0} 3 | Set Ī.sub.r.sup.viol = {i ε Ī.sub.r | λ.sub.r,i > 0} 4 end