Quadratic program solver for MPC using variable ordering

Abstract

A system and approach for storing factors in a quadratic programming solver of an embedded model predictive control platform. The solver may be connected to an optimization model which may be connected to a factorization module. The factorization module may incorporate a memory containing saved factors that may be connected to a factor search mechanism to find a nearest stored factor in the memory. A factor update unit may be connected to the factor search mechanism to obtain the nearest stored factor to perform a factor update. The factorization module may provide variable ordering to reduce a number of factors that need to be stored to permit the factors to be updated at zero floating point operations per unit of time.

Description

BACKGROUND

The present disclosure pertains to control systems and particularly to model predictive control. More particularly, the disclosure pertains to quadratic programming to solvers.

SUMMARY

The disclosure reveals a system and approach for storing factors in a quadratic programming solver of an embedded model predictive control platform. The solver may be connected to an optimization module which may be connected to a factorization module. The factorization module may incorporate a memory containing saved factors that may be connected to a factor search mechanism to find a nearest stored factor in the memory. A factor update unit may be connected to the factor search mechanism to obtain the nearest stored factor to perform a factor update. The factorization module may provide variable ordering to reduce a number of factors that need to be stored to permit the factors to be updated at zero floating point operations per unit of time.

BRIEF DESCRIPTION OF THE DRAWING

FIG. 1 is a diagram of a formula that illustrates box constrained quadratic programming;

FIG. 2 is a diagram revealing a process incorporating variable ordering used when storing factors;

FIG. 3 is a diagram of a module in a context of an application which is integrated with an operating system and an interface;

FIG. 4 is a diagram of an embedded platform that may contain an MPC with a semi-explicit quadratic programming solver having a factorization unit with saved factors, factor search and factor update modules;

FIGS. 5 and 6 are diagrams that reveal an off-line portion and an on-line portion, respectively, of the present approach;

FIGS. 7 and 8 are diagrams that relate to variable ordering finding;

FIG. 9 is a diagram of a table pertinent to selection of factors for storage;

FIG. 10 is a diagram of an illustrative example of a factor storing scheme; and

FIG. 11 is a diagram of a table representing an illustrative example of an ordering to scheme.

DESCRIPTION

The present system and approach may incorporate one or more processors, computers, controllers, user interfaces, wireless and/or wire connections, and/or the like, in an implementation described and/or shown herein.

This description may provide one or more illustrative and specific examples or ways of implementing the present system and approach. There may be numerous other examples or ways of implementing the system and approach.

The present approach may be used for solving QP for embedded MPC application using factor updating. Advanced control problems may usually be formulated as mathematical optimization. As an example, one may mention the model predictive control (MPC) which may often be formulated as a parametric quadratic programming (QP) issue. In the MPC, it may be necessary to solve the parametric QP problem on-line and therefore there may be a need for reliable and fast numerical solvers for QP.

For standard MPC applications in process control, sampling periods may often be in the order of seconds or minutes. Such sampling periods may be long enough to solve a needed QP problem by using a standard powerful PC. MPC may often be used in other than process control applications, such as automotive and aircraft control systems. In the latter applications, sampling frequencies may be higher and the computational resources appear to be limited (CPU, memory). Therefore, if one would like to utilize the MPC control approach under such conditions, there may be a need for fast and tailored QP solvers for embedded applications with limited CPU power and memory.

The present approach appears suitable for solving small QP optimization problems in MPC applications with simple (or box) constraints. Simple constraints for optimization variables may be defined as their lower and upper bounds, “LB<=X<=UB” where X is an optimization variable, LB and UB are the lower and upper bounds, respectively. It may be used for embedded systems (e.g., automotive, aircraft, and so on).

There may be efficient solvers for solving the QP problems with simple constraints. Two basic algorithm classes may be active set solvers (ASS) and interior point solvers to (IPS). The IPS approaches may be suitable for large-scale problems while the ASS approaches may be faster for small to medium size issues. The QP problems in MPC control may be classified as small to medium size problems and therefore ASS approaches appear to be useable. The solvers may be iterative and work with a set of active constraints in each of the iterations. The set of active constraints may be updated during the iterations. During each iteration of an ASS based solver, most of the computational time may be spent on a computation of the Newton step which is in a complexity of N{circumflex over ( )}3 flops (floating point operations per second).

A basic approach may enable one to add or remove a single constraint in a working set which may not necessarily be efficient in a case where there should be large change in the active set to find the optimum.

There may be a class of ASS algorithms for solving the QP problems with simple constraints based only on a gradient projection. The gradient projection based approaches may enable one to add or remove multiple constraints from a working set in each iteration, which can enable quicker identification of a set of active constraints in the optimum. However, these approaches may use a steepest descent direction which might be inefficient and thus the convergence may be too slow, and therefore many iterations may be needed. On the other hand, each iteration may be very inexpensive since only the gradient has to be evaluated in a complexity of N{circumflex over ( )}2 flops.

Finally, there may also exist a class of algorithms which combines both previously mentioned items. The algorithms may use a gradient projection approach for quick identification of the optimal working set, and use the ASS's Newton step computation to improve the convergency, so as to decrease a number of needed iterations. However, these algorithms may involve a computation of the Newton step with a complexity of N{circumflex over ( )}3 flops in each iteration which can lead to slow computation.

The present approach may enable a decrease of the computation time at each iteration by precomputing some part of the Newton step process in combination with a gradient/Newton projection algorithm or an ASS approach or algorithm. In the ASS based algorithm, a KKT (Karush-Kuhn-Tucker) system for currently active constraints may have to be solved at each iteration by using some factorization (LDL, Cholesky, or the like). Thus, it may be possible to precompute and store all or only some factors (partially or whole) of all potential KKT matrices. During a solution, these (partial) factors may be loaded from memory, the factorization process may be finished and then used in the computation of the Newton step in a standard way.

An issue with this approach is that the number of all possible factors of all possible KKT matrices may grow rapidly with the number of QP variables (as with 2{circumflex over ( )}N). This growth may prevent the use of such algorithm for a relatively large N.

The present approach may be similar to a standard way to precompute some factors and store them and use them in the same way in the solution.

But contrary to the standard way, where the factor is computed from scratch if a combination of active constraints corresponding to an unsaved factor is encountered during the on-line computation; in the present approach, the “nearest factor” may be found in the memory and, after that, updated to obtain a factor which is afterwards used in the computation. The “nearest factor” in this sense may be the factor for which the size of the following set is smallest: {N\S}+{S\N}, where N is a set of indices of inactive constraints for a wanted factor, and S is a set of indices of inactive constraints corresponding to a stored factor. The update process may involve additional computation with complexity of N{circumflex over ( )}2 flops which is not necessarily limiting. Also, this approach may lead to the fact that radically less of the factors have to be stored and thus improve the memory limitation of the standard way.

The present approach may be implemented as a software module and be used as a core optimization solver for MPC control for embedded applications (e.g., automotive or aircraft control systems). The approach may be suitable, namely in combination with a gradient projection based QP solver.

As it was said already, one approach for solving a QP problem originated from MPC may incorporate a use of online solvers based on an active set strategy where for each iteration, the so-called Newton process may be computed for a currently estimated active set of working constraints. This may be done via solving a set of linear equations (KKT conditions) usually by a factorization approach which appears as the most computational complex part of the algorithm.

It may be possible to pre-compute and save virtually all factors of possible to combinations of a working set, and then, during the online phase, only load them from a memory to virtually completely reduce the computation burden associated with factorization. On the other hand, such approach may lead to a significant memory need which grows exponentially with a number of optimization variables.

An approach may be to exploit the fact that for box constrained OP issues, the reduced Hessian (THZ) which needs to be factored, can be obtained only be removing the rows and columns of Hessian H which corresponds to the active constraints.

Then when having the corresponding factor of such reduced Hessian with some active constraints, the factor of a new reduced Hessian for the same active constraints plus an additional one may be easily computed only by cutting off the row and column of the original reduced Hessian, i.e., by an updating process. This updating process appears much less computationally expensive (i.e., O(n{circumflex over ( )}2) instead of O(n{circumflex over ( )}3) FLOPS for a computation from scratch).

However, if a special variable ordering is used, then an updating process may be reduced to 0 FLOPS by only removing the last row and column of the factors.

Thus, a goal may be to find a variable ordering in such a way that only a few factors are stored and loaded from memory when needed and the rest of them are created from them by a “fast updating process” with 0 FLOPS complexity. This may limit the computational complexity of the factorization process and also greatly limit the memory need as compared to a case when all of the factors are stored.

A special variable ordering scheme may divide the factors according their size (number of rows) into groups and try to maximize an overlay of variable indices corresponding to inactive constraints at the beginning between the groups, i.e., to allow a multiple cut-off of rows and columns with 0 FLOPS update.

A first approach may be to precompute some factors off-line and save them and then in the on-line phase, use a factor update of the “nearest” factor for currently inactive constraints W. A factor update process may have a complexity O(n{circumflex over ( )}2) per each change—adding or removing one or more rows and/or columns.

For each saved factor, it may be saved information to which a combination of inactive constraints corresponds (such as S).

A nearest factor may be defined as the one with smallest needed changes in the N that is needed to obtain S, i.e., how many rows have to be removed from/added to the stored factor by a factor update process to obtain a factor for currently inactive constraints (those in the set N).

A second approach may be to use an ordering of variables for each combination of indices in N, such that the update process can be done with zero computational complexity first by removing rows of the saved factor.

Ak, Wk may represent a working set of constraint indices indicating which components of optimization vector are on the limits at a k-th iteration of an algorithm. Nk may represent a set of indices of constraints indicating which components of an optimization vector are free, that is not on a limit at the k-th iteration of the algorithm.

FIG. 1 is a diagram that may regard an issue to be resolved, which can be a box constrained quadratic programming as shown by formula 21 with its constraints, having a positive definite Hessian. Such issue may arise in various situations such as an MPC.

FIG. 2 is a diagram of the, present approach. A feature may incorporate is variable ordering to further reduce the factors that need to be stored allowing them to be updated in 0 FLOPs cost. At symbol 81, a set of factors of a KKT system corresponding to an only set of combinations of currently active constraints, may be precomputed and saved. Other factors may be deduced from the saved set by a factor update scheme according to symbol 82. Variable ordering (i.e., permutation matrix) may be used when storing the factors at symbol 83. Deduction of other factors in zero FLOPs just by row and column removal may be enabled according to symbol 84. A factor update scheme may be used with a non-zero cost in FLOPS to further reduce memory consumption at symbol 85. According to symbol 86, the present approach may be suitable for an active set method or a combination of a Newton step/gradient projection algorithm. The approach may also be suitable for fast embedded model predictive control as noted at symbol 87. The items of symbols 81-87 may or may not be noted in numerical order of the symbols.

FIG. 3 is a diagram of a module 13 in a context of an application 14 which is integrated with an operating system 15 and an interface 16. A user or operator 17 may interact with interface 16.

FIG. 4 is a diagram of an embedded platform 21 that may contain an MPC 16. The present approach may be placed within a higher picture relative to the surrounding systems and main idea can be captured there. The present approach may “live inside” some semi-explicit QP solver and consist of three main parts such as a memory with saved factors and algorithm with a search for a factor and a factor update algorithm.

Controller 22 may incorporate a state observer 23 connected to a semi-explicit QP solver 24. Solver 24 may have a QP optimization module 25 that encompasses a KKT system and matrix 27. Solver 24 may also have a factorization unit 26 connected to QP optimization module 25. Factorization unit 26 may have a memory 28 that contains a saved factors module 31. Saved factors module 31 may be connected to a factor search module 30 that is connected to a factor update module 29. Factor search module 30 and Factor update module may be outside of memory 28 but within factorization unit 26. Modules 31, 30 and 29 are a feature of the present approach.

An output of embedded platform 22 may be connected to a unit 32. Unit 32 may have one or more actuators 33 connected to MPC controller 22 and to a physical plant 34 in unit 22, for instance, an engine and/or an after treatment system. Unit 32 may also have one or more sensors 35 connected to physical plant 34 and to MPC controller 22.

FIG. 5 and FIG. 6 are diagrams that reveal an off-line portion 41 and an on-line portion 42, respectively, of an algorithm for the present platform 21. Off-line portion 41 may incorporate an action to find an ordering of variables at each combination of a working set as indicated in symbol 43. The content of symbol 43 may be a feature of the present approach. To select and compute a set of factors for a given MPC QP issue as indicated in symbol 44 may follow the action of symbol 43. Some set of factors may be computed here and only part of them needs to be stored as noted in symbol 45. Thus, symbol 45, which may follow symbol 44, may indicate an action to store a selected sub-set of pre-computed factors to a memory for on-line portion 42.

On-line portion 42 in a diagram of FIG. 5 may be of a semi-explicit QP solver like that of solver 24 in FIG. 3. Portion 42 may begin at symbol 47 followed by an initialization at symbol 48. Following symbol 48, a question of whether to terminate portion 42 may be asked at symbol 49. If an answer is yes, then portion 42 may stop at symbol 51. If the answer is no, then at symbol 52, a KKT matrix may be built based on a current set of active constraints; a nearest stored factor may be found in the memory and a factor update may be performed. The finding the nearest stored factor in the memory and to performing the factor update are feature of the present approach. After symbol 52, a KKT system may be solved and a solution may be updated at symbol 53. Then a working set of active constraints may be updated at symbol 54 and activity of portion 42 may continue on to symbol 49 where the question of termination may be asked.

The present approach in the diagrams of FIG. 4 through FIG. 6 may also be used for a combined Newton/gradient projection algorithm.

FIG. 7 and FIG. 8 are diagrams that relate to variable ordering finding. They may indicate what is variable ordering and how it effects a factor structure. At each iteration, a KKT system may be solved, an equation 57 and an equation 58 may be solved for the KKT system. The terms relative to the equations may be as the following. Z is base of null space of Jacobian of currently active constraints corresponding to the set of indices W, g is a current gradient, pz, is a Newton step in the free components, P is permutation matrix (variable ordering), and p is a Newton step.

P (variable ordering) may be found such that an overlay of Z′HZ is as large as possible for all combinations of active constraints such that only a subset of all possible factors of Z′HZ needs to be solved. This ordering could be found by brute force optimization or by some heuristics.

A diagram in FIG. 8 may reveal an overlay 61 of a factor 1 over a factor 2. Hence Factor 1 can be deduced from Factor 1 by simply removing last row of Factor 2. Another diagram of FIG. 8 may reveal a table 62 an example of possible ordering of variables within factors 1 and 2 such that for inactive set N=I/W.

FIG. 9 is a diagram of a table 64 relative to selection of factors to store. Table 64 indicates an example of variable ordering for problem with three variables. The arrows may stand for “deduction from”. The circled numbers may be saved. Rest of the numbers may be created by cutting them off (i.e., no FLOPS).

Only large factors (corresponding to large number of inactive constraints) may be stored. Variable ordering may be used inside the factors. Zero flops deduction of other factors may be enabled. For illustrative examples, a factor with an inactive set {1 2} may be deduced from factor {1 2 3} only by removal of last row, and a factor with an inactive set {1 3} may have to be solved or stored since there is no factor which enables a 0 FLOPS deduction.

An objective is to speed-up a solution of a KKT system solution. There may be a use of a factor update of a reduced Hessian and variable ordering. By using “clever” variable ordering, then only a subset of factors for each combination of active set may have to be solved. Rest of factors may be computed with 0 flops (floating point operations) by removal of rows.

FIG. 10 is a diagram that illustrates an example of a factor storing scheme. Triangle 71 represents a saved factor for inactive constraints {1,2,3}. Triangle 72 represents a computed factor for inactive constraints {1,2}. Just the last row is deleted as indicated by the shaded area. Triangle 73 represents a computerized factor for inactive constraints {1}. Just the last two rows and columns are deleted. A sequence may be from triangle 71, though triangle 72 to triangle 73. The diagram reveals a basis for a 0 FLOP fast factor update.

FIG. 11 is a diagram of a table 75 representing an example of an ordering scheme for QP problem with six optimization variables (n=6). Indices may indicate which constraints are inactive. An order may then indicate an ordering of corresponding data in a factor. Arrows 76 may indicate the transitions from one index to a lower index. Circles 77 may indicate the numerals that are saved. Numerals may be created by a cutting off (i.e., no FLOPS). Some numerals (not shown) may be computed by an augmented update.

A nearest factor search may be conducted. In an on-line phase the “nearest” stored factor to the currently active constraints W may have to be found. For each saved factor, it is saved information to which combination of active constraints it corresponds (i.e., S). The nearest factor may be defined as the one with smallest needed changes in the Wk is needed to obtain S, i.e., how many rows has to be removed from/added to stored factor by factor update process to obtain factor for currently active constraints (those in the set W).

To recap, a system for quadratic programming may incorporate an embedded platform comprising a model predictive control (MPC) controller connected to a physical subsystem. The MPC controller may incorporate a state observer and a semi-explicit quadratic programming (QP) solver connected to the state observer. The semi-explicit QP solver may incorporate an optimization module and a factorization module connected to the QP optimization module. The factorization module may incorporate a memory having a saved factors unit, a factor search mechanism connected to the saved factors unit, and a factor update mechanism connected to the factor search mechanism.

The factorization module may provide variable ordering to reduce factors which need to be stored to allow them to be updated at a zero floating-point operations per unit time (FLOPs) cost.

The physical subsystem may incorporate a physical plant, one or more sensors situated at the physical plant and connected to the MPC controller, and one or more actuators attached to the physical plant and connected to the MPC controller. The physical plant may be an internal combustion engine or after-treatment device.

The optimization module may incorporate a Karush-Kuhn-Tucker (KKT) subsystem having a KKT matrix.

An approach for storing factors, may incorporate precomputing a set of factors of a matrix corresponding to an only set of combinations of currently active constraints, saving the set of factors, deducing other factors from one or more factors of the set of factors with a factor update scheme, storing the factors with variable ordering, enabling a deduction of the other factors in zero floating-point operations per unit time (FLOPs) by just a row and/or column removal, and using a factor update with non-zero cost in FLOPs to further reduce memory consumption.

The approach may further incorporate computing an active set based approach using pre-computed factors such as a Newton step/gradient projection approach.

The factor update may be used in embedded model predictive control.

The variable ordering may be accomplished with a permutation matrix.

The matrix may be of a Karush-Kuhn-Tucker (KKT) system.

A process for solving a quadratic programming issue, may incorporate an off-line portion, and an on-line portion. The off-line portion may incorporate finding an ordering of variables at each combination of a working set, selecting and computing a set of factors for a model predictive control (MPC) quadratic programming (QP) issue, and storing a selected subset of the set of factors. The on-line portion may incorporate initializing an on-line process, building a matrix on a current set of constraints, finding a nearest stored factor from the memory, performing a factor update with the nearest stored factor, solving the matrix, and updating a solution for a QP issue.

The process may further incorporate updating a working set of active constraints.

The matrix may be a Karush-Kuhn-Tucker (KKT) matrix.

The process may further incorporate reiterating the on-line portion. The process may further incorporate terminating the on-line portion.

U.S. Pat. No. 8,504,175, issued Aug. 6, 2013, and entitled “Using Model Predictive Control to Optimize Variable Trajectories and System Control”, is hereby incorporated by reference. U.S. Pat. No. 8,924,331, issued Dec. 30, 2014, and entitled “System and Method for Solving Quadratic Programming Problems with Bound Constraints Utilizing a Semi-Explicit Quadratic Programming Solver”, is hereby incorporated by reference.

Any publication or patent document noted herein is hereby incorporated by reference to the same extent as if each individual publication or patent document was specifically and individually indicated to be incorporated by reference.

In the present specification, some of the matter may be of a hypothetical or prophetic nature although stated in another manner or tense.

Although the present system and/or approach has been described with respect to at least one illustrative example, many variations and modifications will become apparent to those skilled in the art upon reading the specification. It is therefore the intention that the appended claims be interpreted as broadly as possible in view of the related art to include all such variations and modifications.

Claims

1. A method for storing factors, comprising: precomputing a first set of factors of a matrix corresponding to an only set of combinations of currently active constraints;saving the first set of factors;deducing second factors from one or more factors of the first set of factors with a factor update scheme;storing the first set of factors and the second factors with variable ordering;enabling a deduction of third factors in zero floating-point operations per unit time (FLOPs) by just a row and/or column removal; andusing a factor update with non-zero cost in FLOPs to further reduce memory consumption.

2. The method of claim 1, further comprising computing an active set based method using a Newton step/gradient projection approach.

3. The method of claim 1, wherein the factor update is used in embedded model predictive control.

4. The method of claim 1, wherein the variable ordering is accomplished with a permutation matrix.

5. The method of claim 4, wherein the matrix is of a Karush-Kuhn-Tucker (KKT) system.