Real-time quadratic programming for control of dynamical systems

Abstract

Real-time control of a dynamical system is provided by determining control variables that get as close as possible to producing a desired response. Additional consideration of physical limitations leads to a convex Quadratic Program with inequality constraints that needs to be solved in real-time. A new active set algorithm is described to solve the convex Quadratic Program efficiently that meets real-time requirements. Based on the key observation that the physical limitations of the system translate to optimal active sets that remain relatively unchanged over time (even though the actual optimal controls may be varying), starting guesses for the active set obtained from the final iterate in the previous time period greatly reduces the number of iterations and hence allows the Quadratic Programs to be solved to convergence in real-time.

Description

BACKGROUND OF THE INVENTION

[0002] The present invention relates generally to control systems and more particularly to dynamic inversion control systems using quadratic programming.

[0003] Time varying systems can be modeled by a set of linearized algebraic equations that relate the outputs and performance indices of the system to the inputs at every instant in time (in a nonlinear dynamical system, these equations change at every time instant). These inputs are formally known as control variables, because the values of these variables are under the control of the operator. The aim is to find the best sequence of control variables in time that maximize the performance indices of the system. In addition, the physical constraints of the system must be honored in order to prevent critical system failure. The problem of determining the best sequence of control variables requires solving a sequence of optimization problems in time, i.e., a new optimization problem at every instant in time. Each such optimization problem must be solved almost instantaneously. For the Dynamic Inversion application in systems with fast dynamics, this allows only a few milliseconds. Hence the requirement that the optimization technique come up with a solution swiftly is key in these applications.

[0004] Such a sequence of problems has directly or indirectly been posed for controlling dynamical systems in the past. However, the techniques for solving them have either not operated in real-time, because the applications were not real-time critical or operated in real-time, but the real-time interval was far more generous (perhaps of the order of hours, not milliseconds), as in control of downstream petroleum refining operations and in portfolio optimization. Other techniques operated in real-time, but the technique was such that no rigorous claim could be made about the optimality of the solution, and no attention was paid to not violating the physical constraints of the system. Examples of this technique include certain flight control algorithms.

SUMMARY OF THE INVENTION

[0005] This invention provides an optimization algorithm, which is implemented on an onboard computer for the control of a dynamical system, while honoring the physical constraints of the system to avoid critical failure or breakdown. It is particularly meant for control of systems with relatively quick dynamics, which need quick solution for control variables. The algorithm has special features that make it particularly useful for realtime control applications, as follows.

[0006] The inventive algorithm provides onboard optimization where the solution is restricted to a convex form, in general, and sometimes to quadratic programming, in order to ensure that a unique optimal solution exists, there are no local minima or maxima, and that the solution can be found with reasonable computation.

[0007] The algorithm also uses the active set method for solving realtime convex optimization problems, where the active set of a prior solution is “fedback” to initialize the search for the best active set in the current problem. This envisions application in a control system where the control solution is periodically updated at a high enough frequency that there is only little change from one update to the next.

[0008] The invention also provides an active set method of solution where the iterative search for the best active set is formulated so that every iteration is an improvement over the previous one. This implies that the algorithm can be interrupted, if necessary, or that a given iteration might fail due to unforeseen numerical problems and a good solution can still be provided from a prior iteration.

[0009] The algorithm also uses the active set method, with active set feedback, where the number of iterations in the search for best active set is fixed in order to provide predictable computation time, as required by realtime control software. The linear set of equations are solved using a combination of Cholesky and QR factorizations, in order to provide good numerical properties despite possible degeneracy in the constraints.

[0010] The central challenge is to solve the Quadratic Programming problem within the given real-time interval in which the system is considered quasi-static. In some fast system applications this time interval is only {fraction (1/40)}th of a second, with only a few milliseconds allocated for solving the QP in the Dynamic Inversion. All existing (non-real-time) algorithms for solving QPs are iterative, in that every successive iteration leads to a (not always) better guess for the actual solution. The framework of Accelerated QP is offered by a Quadratic Programming technique based on guessing the set of binding constraints at optimality (the active set), as opposed to an interior point method, for reasons that will become clear. In this Accelerated QP algorithm, every iteration does improve the objective function monotonically. Each of these iterations involves linear algebra computation and takes a certain number of flops (floating point operations), which in turn consumes CPU time. In order to always solve a QP within a fixed (real) time interval, it would thus be necessary to solve the QP within a fixed number of iterations. While there is a theoretical maximum number of iterations within which a solution will be found using an active set method, it is not possible to always execute those iterations within the given time interval using the existing controller hardware.

[0011] Thus Accelerated QP focuses on solving each QP in as few iterations as possible, fixing the number of iterations or computation time, and making each QP iteration computationally so efficient so that each iteration takes up minimal CPU time.

[0012] Key ideas behind achieving minimal iterations:

[0013] At each iteration, the active set QP algorithm guesses which constraints are binding at the solution. This guess is based on the signs of the Lagrange multipliers of the constraints estimated to be binding at that iteration (called the working set of constraints). A negative sign Lagrange multiplier indicates that it is possible to move in the interior of the corresponding constraint and improve the value of the cost function. While traditional QP algorithms drop only the constraint with the most negative Lagrange multiplier from the working set, Accelerated QP drops all the constraints with negative multipliers from the working set. This is because the traditional method may end up dropping the constraints one by one, leading to more iterations. Dropping all the constraints with negative multipliers allows the next QP iteration freedom to move in a larger subspace, and thus allows more rapid convergence to the solution in fewer iterations. This is of key importance, because in a real-time setting where the operations must be completed within a certain time interval, the QP is stopped after a fixed number of iterations.

[0014] The other important idea in minimizing the number of QP iterations required for convergence is that of a ‘hot start.’ A sequence of Quadratic Programs is solved in time, as required by the evolving dynamical system. Since the constraints in the QP are related to the physical state of the system, the constraints that are binding at the solution directly correspond to the physical factors limiting the system. Since these limiting factors change only little, or not at all, from one time interval to the next, the set of binding constraints from the QP in the current time step is an excellent, and most often exact, predictor of the set of binding constraints at optimality for the QP in the following time step (even though the parameters defining the QP are different). This key observation allows us to solve all the QPs to optimality in only a few iterations, which is 3-5 iterations in the worst case for a Matlab simulation for Dynamic Inversion.

[0015] If the short real time interval does not allow the solution of a QP to reach optimality, it still exits with the best solution obtained thus far. Now, since the QP that is not solved to optimality is almost similar to the QP in the following time step, the working set for the non-optimal solution is quite close to the optimal binding set for the next QP. Thus, even though the solution to a QP in a time step may not be optimal, eventually a time step is reached when the QP is solved to optimality, because the starting guess is so good. In this manner, the ‘hot start’ scheme allows the algorithm to ‘recover from non-optimality’ caused by real-time limitations in a matter of a few time steps. In computer simulations for Dynamic Inversion, optimality is seen to be achieved in just one extra time step.

[0016] The active set method is naturally capable of exploiting the ‘hot start’ information, which is why it forms the foundation of Accelerated QP, as opposed to an Interior Point Method.

[0017] Key ideas behind speeding up each QP iteration:

[0018] The ideas behind speeding up each iteration of the active set search rely on a combination of efficient linear algebra and performance programming methods, as is described in detail in an upcoming section

BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Other advantages of the present invention will be readily appreciated as the same becomes better understood by reference to the following detailed description when considered in connection with the accompanying drawings wherein:

[0020] FIG. 1 illustrates the multivariable control system of the present invention.

[0021] FIG. 2 shows graphs of the response of the system of FIG. 1 when no limits are reached.

[0022] FIG. 3 shows graphs of the response of the system when limits are reached.

[0023] FIG. 4 shows the computational burden of the QP solver of FIG. 1.

DESCRIPTION OF THE PREFERRED EMBODIMENT

[0024] FIG. 1 illustrates a control system 20 implementing the dynamic inversion via real-time quadratic programming algorithm of the present invention. The control system 20 is briefly summarized here and described in more detail in co-pending application U.S. Ser. No. ______, entitled “Constrained Dynamic Inversion Control Algorithm,” filed on the same day and by the same inventors as this patent application, and which is hereby incorporated by reference in its entirety. The constrained dynamic inversion algorithm is meant to be implemented on a discrete time sampled-data computer 22 in order to control a dynamic system. The computer 22 includes a CPU, memory 24 and other known hardware and is programmed to perform as described herein to produce actuator commands, uc. The actuator commands, uc, control a plurality of actuators 26, which take the electronic commands, uc, and turn them into physical forces that adjust a plurality of effectors 28, each of the plurality of actuators 26 being associated with one of the plurality of effectors 28. The effectors 28 change the dynamics of the plant 30. A plurality of sensors 32 provide feedback y to the computer 22. The computer 22 determines the next set of actuator commands, uc, based upon the information from the sensors 32 and based upon system output commands 34, which may come from a human operator or a higher level control system.

[0025] The computer 22 periodically samples commands, yc, and sensor feedback, y, in order to make actuator commands, uc. y(n) and yc(n) refer to values of y and yc at the nth sampling, which occurs at time n. The actuators 26 interpret uc(n) as a constant command held from time n to time n+1. y, yc, and u can all be vectors to represent multiple degrees of freedom. As used herein, the current time is time n and the time of the next controller sample is time n+1. The control algorithm is implemented via controller software. It determines the actuator commands, u, that make the system outputs, y, match the commanded values yc. The commands may come from a human operator or a higher level control system. This is a well known form of a digital computer based control system. The substance of this invention is a novel control algorithm stored in memory 24 and executed by computer 22.

[0026] The control logic (software) architecture separates control logic design and implementation into Dynamic Feed-Forward (DFF) 40, Dynamic Model Follower (DMF) 42 and Dynamic Inversion (DI) 44 modules. The DFF 40 module produces a component of a command of how the system output should change by the next controller sample period. Its design establishes the nominal command response dynamics. The DI 44 module determines effector 28 commands that best cause the system outputs to change as desired. DI 44 thus “inverts” the system dynamics by making them appear as simple, uncoupled digital integrators. Even when the system dynamics vary with configuration or operating point, the dynamics with DI 44 appear the same throughout the system's operating envelope and the same from one type of system to another. The DMF 42 module compares desired and actual system responses to adjust the output change command as needed to attenuate the effects of modeling errors and disturbances. DMF 42 establishes how quickly and completely the effects of model errors and disturbances are attenuated.

[0027] The DI 44 module can be further subdivided into four modules as shown in FIG. 1, for a total six modules as follows:

[0028] A state variable model 48 produces a linearized model for predicting values of system variables, y(n+1), one control sample ahead (a.k.a. the next step) as a function of the actuator commands, uc(n). This prediction also uses a feedback measure of the current output, y(n).

[0029] A performance index block 50 uses this model to formulate the determination of actuator commands, uc(n), as a quadratic optimization problem with linear constraints. One input to this block is the command v(n) which is the amount the system output y should change by the time n+1. The performance index contains a quadratic measure of how well this will be achieved, as predicted by the model.

[0030] A constraint block 52 transforms limits and constraints on the effectors and on various variables of system response into the form of linear inequality equations with uc(n) as the independent variables.

[0031] A quadratic programming block 54 solves the special type of optimization problem poised by the performance index and constraint equations. The quadratic programming block 54 is the focus of the present invention.

[0032] Optimization-based controllers offer several advantages. First, they offer a general, reusable method for enforcing constraints in a multivariable system. The general mathematics of constrained optimization will cause the system behavior to change as constraints are met. Secondly, they are flexible. The model, constraints, and performance index can vary between controller sample times. Control actuator failures and other system failures can be reflected in the controller as changes to the model and the imposed constraints. Once these change, the optimization-based controller will attempt to do the best job it can with the functional controls. Finally, optimization based controllers can operate at a higher level of autonomy. They can receive commands in the form of numerous and possibly conflicting goals with various and shifting priorities. The optimization enables a tradeoff among the goals. It is important to note that providing these advantages does not require the precise optimal solution per se. Instead, optimization is being used as a mathematical framework for designing control software at a higher level of abstract so it is more flexible and reusable.

[0033] Controllers operate automatically in hard real time. Because of this, associated optimization software must be highly reliable and must complete its computations within a fixed time interval. For control theoretic reasons, this interval is generally a twentieth to a tenth of the fastest time constant relevant to system operation. For example, a high performance aircraft needs to substantially attenuate the effect of a sudden gust within 0.3 seconds. This implies the controller needs to complete its computations in approximately 25 msecs. Therefore, the optimization software needs to be computationally efficient to be practically implementable.

[0034] The high level strategy for real time constrained optimization is to make the approximations necessary to pose the optimization as convex quadratic programming (QP) problems. Every convex QP problem has a unique global optimum. Algorithms for solving a QP reliable search procedures, involving predictable and reliable numerical properties. A QP problem can be formed by measuring performance with quadratic indices, expressing limits as locally linear inequality constraints, and modeling the system with locally linear, possibly parameter varying, difference equations. These are mild approximations in most applications.

[0035] The description below focuses on the formulation and solution of an instance of such a QP at every point in time, given a fixed budget for floating point operations. Clearly, it is difficult to promise that the algorithm will arrive at a solution within a fixed number of floating point operations (unless the number of floating point operations is so large that it allows a number of iterations equal to the number of combinatorial possibilities for working sets). Thus, special emphasis is given on arriving at a very good solution in the allowed time or within the floating point operation (flop) budget.

[0036] Problem Description

[0037] Consider a dynamical system with control variables u(t), u: →n, and state variables ξ(t), ξ:→p. Consider a linear model for the dynamics of the system: 1$\frac{d\xi }{dt}=A_1\xi +B_{1}u$$\frac{du}{dt}=B_{2}u+M{u}_c.$

[0038] Here uc represents the control command issued to the controllers. The responses of the system, y(t), y:→p are defined as linear combinations of the state and control variables:

y (t)=Cξ(t)+Du(t)

[0039] The objective is to find uc(t) so that dy/dt matches a desired vector v as best as possible. This can be posed as 2$\min \&LeftDoubleBracketingBar;\frac{dy}{dt}-v\&RightDoubleBracketingBar;w$

[0040] where ∥.∥w represents a weighted l2 norm, often referred to as an elliptic norm, where W is a diagonal matrix of positive weights. The above linear system can be discretized (with Δt=1) as

ξt−ξt−1=A1ξt−1+B1ut−1  (1)

y t −y t−1 =C (ξt−ξt−1)+D(ut−ut−1)  (2)

[0041] Under the assumption that controller dynamics is “fast”, i.e. B2=−I and M=I, so that uc=ut. This is done mainly for illustrative purposes to simplify the algebra, and is not fundamentally any more difficult to handle. Substituting 1 in 2 yields

y t −y t−1 =C (A1ξt−1+B1ut−1)+D(ut−ut−1)  (3)

[0042] At time point t, the aim is then to minimize

min∥(yt−yt−1)−v∥w

[0043] or, substituting 3, 3$\begin{array}{cc}\underset{u_t}{\min }\&LeftDoubleBracketingBar;{\mathrm{Du}}_t+{\mathrm{CA}}_1{\xi }_{t-1}+\left({\mathrm{CB}}_1-D\right)u_{t-1}-v\&RightDoubleBracketingBar;w& \left(4\right)\end{array}$

[0044] Note that in the above objective function at time t, the quantities from time t−1 are already known, and ut is the variable the objective is optimized over. This is in contrast with Model Predictive Control where a time window of length greater than one unit of time is used.

[0045] The objective function 4 is a convex quadratic in ut. In addition, there are bound constraints on and possibly other linear constraints involving yt, Δyt, Δut, that prevent the solution of the QP from resulting in actuator commands that may cause the system to fail. These can be rewritten as linear inequality constraints in ut, and added to the mathematical programming formulation.

[0046] It so happens that the system under consideration is often operated under very extreme conditions. Consequently, it is important to take these inequality constraints into account. This is because a simpler solution to the QP without the inequality constraints would, in fact, violate those constraints under extreme conditions and lead the system to critical failure.

[0047] The system design has more actuators, hence control variables, than responses, hence the matrix D is under-determined. The available degrees of freedom are used to attain a secondary objective in meeting desired values, ud, for ut. This secondary object can be expressed as 4$\begin{array}{cc}\underset{u_t}{\min }\&LeftDoubleBracketingBar;u_t-u_d\&RightDoubleBracketingBar;s& \left(5\right)\end{array}$

[0048] The objective in 4 is penalized with the secondary objective above to obtain the combined objective. Squaring and adding the two parts for differentiability, the combined objective function is 5$\begin{array}{cc}\underset{u_t}{\min }{\&LeftDoubleBracketingBar;{\mathrm{Du}}_t+{\mathrm{CA}}_1{\xi }_{t-1}+\left({\mathrm{CB}}_1-D\right)u_{t-1}-v\&RightDoubleBracketingBar;}_W^2+\rho {\&LeftDoubleBracketingBar;u_t-u_d\&RightDoubleBracketingBar;}_S^2& \left(6\right)\end{array}$

[0049] where ρ is a small scalar denoting the relative importance of the second objective, typically 10−3.

[0050] The addition of the secondary objective as a penalty term has a regularizing effect on the first term. It increases the magnitudes of the diagonal elements in the second order term in ut and makes the eigenvalues strictly positive. The QP thus becomes a strictly convex QP.

[0051] Role of secondary objective Having redundant controllers provides additional freedom that can be exploited in order to achieve secondary goals, such as the following:

[0052] Increasing the robustness of the system in the event of an actuator failure.

[0053] Operating the engine at a lower temperature, and hence increasing the life of the engine.

[0054] Maintaining control variable levels close to pre-determined best settings.

[0055] Increasing fuel efficiency of the engine.

[0056] Increasing margins for fast maneuvers. E.g. an actuator that responds faster can be set further away from its bound than a slower actuator. The centered faster actuator can then be used for quick response for a demanding swift maneuver.

[0057] Rename the optimization variables as x instead of ut to avoid cumbersome notation, the reduced problem can be put in the form of the following strictly convex QP with only inequality constraints: 6$\begin{array}{cc}\underset{x}{\min }q\left(x\right)=\frac{1}{2}{x}^T\mathrm{Hx}+c^{T}x& \end{array}$ s.t. Ax≦b  (7)

[0058] The challenge lies in arriving at a good solution within the allowed real time interval in a robust manner.

[0059] Further, H, c, A and b all vary from one time step to the next, but the subscript t is dropped for simplicity of notation (considering only QP at a particular time point here).

[0060] Solving a Convex QP

[0061] The two broad categories of techniques for solving a convex QP differ in how the inequality constraints are handled in the algorithm.

[0062] Active Set Method The older of the two is the active set method, where some of the inequalities are treated as equality constraints at every iteration and the next step is taken such that these equalities are not violated. The set of inequalities which are held as equalities, or are binding or active is referred to as the active set. Feasibility with respect to all the inequality constraints is maintained at every iteration. The Simplex Method for linear programming is a well-known example of an active set method. Theoretically, an active set algorithm could end up visiting all the vertices of the polytope, i.e., end up searching through all possible combinations of the constraints (as demonstrated by the pathological Klee Minty examples for linear programming). However, in practice, convergence is quite rapid.

[0063] Interior-Point Approach The other alternative is an Interior Point method, which can be thought of as applying Newton's method to the perturbed first-order necessary conditions for optimality of the QP. Iterates stay strictly feasible at every iteration (i.e., none of the inequalities are allowed to become binding). For a convex QP, the iterates could stay very close to what is known as the central path in course of convergence. Polynomial time convergence can be proved for Interior Point methods, i.e., they do not demonstrate the worst case exponential time behavior as a standard Simplex method does on the Klee Minty examples. Interior point methods are very well-suited for handling large problems with a large number of inequality constraints.

[0064] Why an active set method is more fitting In this case, some insight is offered by the physical system. The system is very often operated under extreme conditions, and it is not expected that the system to jump from one set of extreme conditions to another set of extreme conditions rapidly. Since the extreme condition in turn govern which inequality constraints are binding at the solution, the set of inequalities that are bonding at the solution to the QP at a particular time step provides a very good, often perfect, estimate of the active set at optimality for the QP in the next time step. This is a good reason to choose the active set method over the interior point method for the problem, since the solution can be used in the current time step to initialize the QP in the next time step.

[0065] An Active Set Method for the Problem

[0066] At the initial time point to, the QP is solved to optimality using an active set method. At any time point t>to, there is a limited budget of flops, which translates into a maximum number of iterations that can at most be executed. The procedure for thus solving the QP is sketched below.

[0067] Find initial feasible starting point. The fact that the active set method produces iterates that are feasible with respect to the constraints relies on having a feasible starting point to begin. For a general QP, one would need to solve a “phase 1” linear program to find a feasible point. Fortunately, the physical system under consideration permits the assumption that x=0 within the feasible set. x=0 corresponds to not issuing any new control commands in the current time point, and thus maintaining the state of the previous time point, which, though perhaps not optimal, should at least satisfy the constraints.

[0068] Subsequent iterations. Form the Working set. While the Active set is the set of inequality constraints that satisfy their right hand sides as equalities at the current iteration, the working set is a subset of the active constraints. This subset does not have to be a strict subset. For the algorithm, up to the first n active constraints are included in the working set, where n is the number of variables. In practice, the number of active constraints is much less that n, and hence all the active constraints are in the working set.

[0069] Given xk, the current kth iterate, determine the next step sk by solving the equality constrained QP 7$\underset{s_k}{\min }q\left(x_k+s_k\right)=\frac{1}{2}{\left(x_k+s_k\right)}^{T}H\left(x_k+s_k\right)+c^T\left(x_k+s_k\right)$ s.t. Ek(xk+sk)=rk

[0070] where Ek is a matrix whose rows correspond to the rows of A containing the constraints in the working set, while the components of rk contain the corresponding subset of right hand side entries from b. In plain words, this above problem looks for the next point xk+sk to be optimal with respect to the quadratic objective, while forcing the constraints in the working set to be held as equalities. However χk+sk may well lie outside the feasible region.

[0071] Since xk is feasible and Ek contains rows that were active in the previous iteration, Ekxk=rk. Substituting this in the above and re-casting the objective leads to the following equality constrained QP for determining the next step. 8$\underset{s_k}{\min }\frac{1}{2}{s}_k^T{\mathrm{Hs}}_k+{\left(c+{\mathrm{Hx}}_k\right)}^T{s}_k$ s.t. Eksk=0  (8)

[0072] The positive definiteness of H allows us to express the unique global minimizer of (8) in closed form, but special attention needs to be paid in the case when Ek is rank-deficient. This will be elaborated in the next section. Assuming H is invertible (which is guaranteed by the addition of secondary objectives with strictly positive weights), the closed form solution to (8) turns out to be

λk*=(EkH−1EkT)−1EkH−1(c+Hxk)  (9)

and

s k *=−H −1 (c+Hxk+EkTλk)  (10)

[0073] where λk denotes the Lagrange multipliers of the constraints Eksk=0.

[0074] The implementation of the solution in this form thus allows us to pre-compute the Cholesky factorization of H (a proxy for its inverse) at the onset of the iterations, and then re-use this already computed factorization in every iteration, thus reducing the overall computational cost per iteration in comparison with solving the original linear system at every iteration. Thus, this procedure exploits the fact that H is invertible in avoiding the solution of a much larger linear system at each iteration.

[0075] Drop constraints corresponding to negative components of the multipliers λk from the working set. A negative value of the multiplier indicates that it is possible to move in the strict interior of the inequality constraint and improve the value of the objective. Hence those inequality constraints cannot be binding at the optimum.

[0076] Given the step sk, determine the maximum value of α∈[0,1] such that χk+αsk is feasible. In other words, the step is scaled back so that the updated point is feasible. This leads to a new set of active constraints.

[0077] Update the current iterate

x k+1 =x k +αs k

[0078] If no constraints were dropped from the working set, i.e. all multipliers were non-negative, and if no new constraint became active, i.e., the step lies within the feasible region and α=1, then the current point xk+1 is optimal.

[0079] Otherwise, return to the initial step to choose the working set.

[0080] While it can be proved that the above iterations would converge to the global optimum in a finite amount of time, a solution must be provided within the allowed time frame. Thus the execution of only a pre-set maximum number of the iterations described above is allowed, depending on how many flops the processor can handle within the allowed real-time interval. Part of the flop budget is also set aside to compute the global optimum of the unconstrained objective scaled back to lie within the feasible region. This point, x*, is computed as follows:

x*=−αH −1 c   (11)

[0081] where α is the largest number in [0, 1] such that x* lies within the feasible region. Such a value of α always exists given the earlier assumption that x=0 is always feasible.

[0082] Now given the above solution, which referred to as the scaled, unconstrained solution, and the best point obtained by the QP iterations, the best solution thus far is returned. The QP algorithm produces feasible iterates that improve the objective at every iteration. A formal proof of this based on convexity is provided in the following sub-section.

[0083] It is possible that the QP may not be solved to optimality in some of the time intervals. However, it is hoped that the final iterations in the current time interval will provide a good lead or a hot start for the QP in the next time interval, and thus convergence will eventually be attained in a subsequent time interval. It does not seem necessary because convergence is attained for every QP in the allowed time interval.

[0084] QP iterates improve the objective function monotonically at every iteration, which can be proved rigorously.

[0085] Handling Degeneracy in the Constraints and Computing the Step and the Multipliers

[0086] Degeneracy is caused by introducing redundant constraints that could alternatively have been expressed as linear combinations of other constraints. Often, such dependent constraints could become active at the same time, and end up in the working set. This results in the loss of row rank in Ek, and needs to be dealt with in the computation. In the simplex method for linear programming, where only one active constraint is added to the working set and only one is dropped in a given iteration, degeneracy could lead to ‘cycling’ amongst some choices of the active sets, and can be rectified with simple rules like Bland's Rule. However, the active set is constructed differently for this QP algorithm, and cycling is not an issue dealt with here. Nevertheless, handling the computation of the multipliers in (9) when Ek loses rank is not necessary, and this is chalked out below.

[0087] The essential idea is to use a QR factorization of EkH−1EkT. This may appear unconventional, since in nonlinear programming, it is more common to perform a QR factorization of ET and then use the factorization to come up with a basis for the null space of E and use it in computing λ and sk. However, this approach requires fewer flops for this particular algorithm.

[0088] Compute wk=c+Hxk

[0089] Use a sequence of back solve and forward solve using the Cholesky factorization (L) of H computed earlier to solve for the columns of the matrix Bk

HB k =E k

[0090] Compute rk=Bkwk

[0091] Compute Mk=BkEkT

[0092] Perform a QR factorization of the m×m matrix Mk, where m is the number of constraints in the working set. So

QR=E k H −1 E k T

[0093] Form QTrk. Since Q is orthogonal, the multipliers can be obtained by solving

Rλ k =−Q T r k

[0094] This is only a forward solve involving O(m2) operations since R is upper triangular. This is where the effects of degeneracy come into play. A degenerate working set can cause the diagonal entries in R corresponding to the dependent constraints to become zero. So whenever a zero pivot is encountered during the forward solve, the corresponding component of λk is set to zero as well. Mathematically, the values of the multipliers can be set to any number. However, they are set to zero to indicate that these multipliers correspond to redundant constraints that do not play a role in the current subproblem.

[0095] Once λk has been determined, the step in x is obtained as

s k =−H −1 (wk+EkTλk)

[0096] Exploiting Performance Programming for Computational Efficiency

[0097] Maximal advantage of computational efficiencies is taken in performing the linear algebra computation in the above formulas in an object-oriented C++ implementation. Along with avoiding common pitfalls like creating copies of objects, frequently called functions such as access operators are automatically inlined, and other features of C++ are used to result in the efficiencies listed below. It is believed that there is no code for Quadratic Programming, available commercially or otherwise, that employs the suite of efficiency-maximizing schemes below:

[0098] While accessing elements of a matrix, it is more efficient to access adjacentd blocks of memory at once. The linear algebra involved in the QP algorithm was carefully analyzed, and it was determined that it was best to access the elements of some matrices row-wise, while for others, the elements had to be accessed column-wise. The pre-defined two dimensional arrays in no programming language would allow for two types of access (e.g., in C, the 2-D arrays are stored row-wise, and in Fortran, columnwise). It is only by means of defining objects in C++ that the Accelerated QP algorithm can make use of matrices that store the elements differently. Based on how the elements are stored, custom routines were also implemented for matrix-vector multiplications that take full advantage of the storage schemes. Other matrix classes such as Upper and Lower Triangular and Symmetric Matrix were also implemented which enable the Accelerated QP algorithm to take full advantage of the storage schemes.

[0099] The QP algorithm uses a matrix whose rows are a sub-set of the rows of the constraint matrix for the original problem. The Accelerated QP implementation in C++ allows for a ‘parasite matrix’ that contains no memory of its own but simply points to the rows of the parent matrix. All linear algebra operations possible on the host matrix are possible on the parasite matrix. This helps to avoid a lot of copying.

[0100] It is necessary to solve a linear system at every QP iteration in order to determine the multipliers for the working constraints, which forms the bulk of the computational cost. For the sake of stability, this is done by means of a QR factorization. The key efficiency enhancement here was in recognizing that the Q matrix need not be formed in order to solve the system. All that is required is multiplication of the transpose of Q by a vector. This can be done by using the Householder transformation vectors for the QR factorization without forming Q at all, thus saving a good amount of flops. Another important thing to realize here is that very often the dimension of the system solved was only 1×1 or 2×2, because the dimension of the system is equal to the number of constraints in the working set. Since the linear system can be solved trivially for the 1×1 case and more directly using a closed form formula based on determinants for the 2×2 case, these two cases are handled separately and thus require many fewer flops than with the QR factorization.

[0101] It is important to emphasize that in developing the above algorithm, primary attention was given to making it stable, dependable and robust. So the Accelerated QP algorithm is not only swift, but also very robust. Contrary to convention, the QR factorization approach is used to compute the Lagrange multipliers in order to handle cases that may be degenerate, i.e., where constraints that become binding simultaneously are linearly dependent. This situation may be caused by redundant constraints that reflect the same physical limitation. The system to which this invention is applied may have such redundant constraints built in for the sake of robustness.

[0102] Further, in the algorithm, some matrix vector multiplications for rows corresponding to constraints in the working set are avoided, because the results are known by virtue of the fact that they are in the working set. This avoids small (but critical) floating point errors that creep up otherwise, while lowering flops as an added bonus.

Claims

1. A method for controlling a multivariable system including the steps of: a) Receiving a plurality of sensor signals indicating current conditions of the system, which are sampled; b) Receiving a plurality of commands; c) Determining a desired dynamic response of the system based upon the commands and the sensor signals; d) in each of a plurality of time steps, formulating a problem of controlling the system to achieve the desired dynamic response as a solution to a quadratic programming problem; e) solving the quadratic programming problem in each time step using an iterative algorithm which searches for an optimal active set; f) in each subsequent time step of the plurality of time steps, initializing a search for the best active set with a final active set of a prior time step of the plurality of time steps.

2. The method of claim 1 wherein the active set comprises a set of constraints that are binding at an optimized solution.

3. The method of claim 2 wherein said step e) further includes the procedure of dropping at least one constraint with a negative Lagrange multipliers.

4. The method of claim 3 wherein said step e) further includes the procedure of dropping all constraint with negative Lagrange multipliers.

5. The method of claim 1 further including the step of using a QR factorization to improve robustness in the event of degeneracy.

6. The method of claim 1 where a total number of iterations of the iterative algorithm is fixed in order to be time predictable.

7. The method of claim 1 further including the step of formulating the quadratic program at every iteration in a manner so that the next iterate is an improvement on a previous iterate.

8. The method of claim 1 further including the step of: precomputing a Cholesky factorization of H; and using the precomputed Cholesky factorization of H in a plurality of iterations of the iterative algorithm of step e).

9. The method of claim 1 furthering including determining the Cholesky factorization via a QR factorization of a generalized square root matrix of H.

10. A multivariable control system comprising: A plurality of sensors indicating current conditions of the system; A plurality of effectors for changing dynamics of the system in response to effector commands; A dynamic feed forward module determining a desired dynamic response; A dynamic model follower comparing desired and actual system responses to adjust an output change command as needed; A dynamic inversion module determining the effector commands in order to change the system output based upon the determination by the dynamic feed forward and dynamic model follower modules by formulating a problem of controlling the system to achieve the desired dynamic response as a solution to a quadratic programming problem in each of a plurality of time steps; and A quadratic programming module finding an optimized solution to the quadratic programming problem in each time step using an iterative algorithm which searches for an optimal active set, the quadratic programming module initializing a search for the best active set in each subsequent time step with a final active set of a prior time step.

11. The system of claim 10 wherein the quadratic programming module guesses the active set based upon Lagrange multipliers of a subset of the constraints.

12. The system of claim 11 wherein the quadratic programming module drops all constraints with negative Lagrange multipliers.

13. The system of claim 10 wherein the active set comprises a set of constraints that are binding at an optimized solution.

14. A computer readable medium storing a computer program, which when executed by a computer performs the steps of: a) Receiving a plurality of sensor signals indicating current conditions of the system, which are sampled; b) Receiving a plurality of commands; c) Determining a desired dynamic response of the system based upon the commands and the sensor signals; d) in each of a plurality of time steps, formulating a problem of controlling the system to achieve the desired dynamic response as a solution to a quadratic programming problem; e) solving the quadratic programming problem in each time step using an iterative algorithm which searches for an optimal active set; f) in each subsequent time step of the plurality of time steps, initializing a search for the best active set with a final active set of a prior time step of the plurality of time steps.

15. The computer readable medium of claim 14 wherein the active set comprises a set of constraints that are binding at the optimized solution.

16. The computer readable medium of claim 15 wherein said step f) further includes the step of guessing the active set based upon Lagrange multipliers of a subset of the constraints.

17. The computer readable medium of claim 16 wherein said step f) further includes the step of dropping all constraints with negative Lagrange multipliers.