Model predictive control with uncertainties

Abstract

An operation of the machine is iteratively controlled with control inputs determined from an optimization of a cost function along a prediction horizon subject to constraints on the control inputs. The optimization is performed according to the model and the cost function includes at least one parameter of the model. During at least some iterations of the control, the parameter of the model is updated resulting in updating the cost function.

Claims

1. A method for controlling an operation of a machine according to a model of the machine having parameters of a physical structure of the machine, comprising: controlling iteratively the operation of the machine with control inputs determined from an optimization of a cost function along a prediction horizon subject to constraints on the control inputs, wherein the optimization is according to the model, wherein the cost function includes at least one parameter of the model, wherein the parameter of the model is uncertain due to an uncertainty of a value of a component of the physical structure of the machine represented by the parameter; and updating the parameter of the model resulting in updating the cost function for at least some iterations of the controlling, wherein steps of the method are performed using a processor, and wherein the controlling changes a state of the machine, wherein the cost function includes a stabilization component for directing a movement of the machine to a target state and a performance component for optimizing the movement of the machine until the target state, and wherein the stabilization component includes an estimate of the parameter of the model, wherein the stabilization component includes a sum of stabilization matrices weighted by the estimate of the parameter of the model, and wherein the parameters of the model are subject to bounded uncertainties defining ranges of admissible values of the parameters of the model, wherein values of the ranges of admissible values form a collection of models of the machine, such that the model of the machine is a convex hull of the collection of the models of the machine, further comprising: determining at least part of the stabilization matrices using a symmetry between values of the extreme models.

2. The method of claim 1, wherein the model of the machine includes $x\left(t+1\right)=\hat{A}x\left(t\right)+\hat{B}u\left(t\right)=\underset{i=1}{\overset{l}{.Math.}}{\left[\left(t\right)\right]}_i\left(A_{i}x\left(t\right)+B_{i}u\left(t\right)\right),$ and the cost function includes ${x\left(N\right)}^{}P\left(\left(N\right)\right)x\left(N\right)+\underset{k=0}{\overset{N-1}{.Math.}}{x\left(k\right)}^{}\mathrm{Qx}\left(k\right)+{u\left(k\right)}^{}\mathrm{Ru}\left(k\right),$ wherein x is the state of the machine, u is the control input, is an uncertain convex combination vector of length l and []i denotes its i.sup.th component, A.sub.i, B.sub.i, i=1, . . . , l are state and input matrices for a collection of models, , {circumflex over (B)} are estimates of the state and input matrices of the machine model, N is a length of the prediction horizon, Q, R are state and input weight matrices, and P() is a sum of stabilization matrices dependent on .

3. The method of claim 1, wherein the parameters of the model are subject to bounded uncertainties defining ranges of admissible values of the parameters of the model, wherein values of the ranges of admissible values form a collection of models of the machine, such that the model of the machine is a convex hull of the collection of models of the machine, further comprising: determining the parameters of the model as weights of the convex hull of the collection of the models of the machine fitting the control inputs and the state of the machine.

4. The method of claim 3, wherein the constraints includes control invariant constraints that ensure satisfaction of specification constraints for the admissible values of the parameters of the model and admissible control inputs to the machine.

5. The method of claim 3, wherein the constraints include a terminal constraint specifying a region for a terminal state at the end of the prediction horizon, further comprising: determining the terminal constraint using specification constraints, the convex hull of the collection of models and the cost function.

6. The method of claim 5, further comprising: determining the region of the terminal constraint using a backward reachable method that iteratively determines, starting from a current set of feasible states of the machine, a backward reachable set of the states within the current set that after a one step evolution of dynamics of the machine is inside the current set until the backward reachable set of states of the machine is enclosed by the current set.

7. The method of claim 4, further comprising: determining a length of the prediction horizon for which there is a solution to the optimization of the cost function for all states satisfying the control invariant constraints; and determining a length of the prediction horizon as largest between a preferred length of the prediction horizon and the determined length of the prediction horizon, wherein the preferred length of the prediction horizon is a predetermined value.

8. The method of claim 7, wherein the length of prediction horizon is equal to the preferred length the terminal set is computed from stabilizing gains and from specification constraints, and the control invariant set is computed from the terminal set and the prediction horizon having the preferred length.

9. The method of claim 1, wherein the prediction horizon includes a sequence of prediction instants, such that the optimization produces control input for each prediction instant forming a sequence of control inputs for the prediction horizon, wherein the parameters of the model for determining the sequence of the control inputs changes between at least some prediction instants.

10. The method of claim 9, further comprising: determining, for each prediction instant the model uses the estimate of the parameter of the model determined a number of time instants before the prediction instant to produce a sequence of prediction models for the prediction horizon with different estimations of the parameter; and jointly determining the sequence of control inputs using the sequence of prediction models.

11. The method of claim 10, wherein the number of time instants is defined by a length of the prediction horizon.

12. The method of claim 1, further comprising: converting the optimization of the cost function into a quadratic program; and solving the quadratic program to determine the control inputs.

13. A system for controlling an operation of a machine according to a model of the machine having parameters of a physical structure of the machine, comprising: a sensor for determining a state of the machine; a processor for updating a value of at least one parameter of the model using the state of the machine, such that the update value of the parameter satisfies a constraint on admissible values of the parameter; and a controller for controlling iteratively the operation of the machine with control inputs determined from an optimization of a cost function along a prediction horizon subject to constraints on the control inputs using the model for prediction along the prediction horizon, wherein the cost function includes at least one parameter of the model, such that the updated value of the parameter updates the model and the cost function, wherein the parameter of the model is uncertain due to an uncertainty of a value of a component of the physical structure of the machine represented by the parameter, wherein the cost function includes a stabilization component for directing a movement of the machine to a target state and a performance component for optimizing the movement of the machine until the target state, and wherein the stabilization component includes a sum of stabilization matrices weighted by estimates of the parameters of the model determined using a symmetry of a convex hull of a collection of models of the machine formed with values of the parameters defined by ranges of admissible values of the parameters of the model.

14. The system of claim 13, wherein parameters of the model are subject to bounded uncertainties defining ranges of admissible values of the parameters of the model, wherein values of the ranges of admissible values form a collection of models of the machine, such that the model of the machine is a convex hull of the collection of models of the machine, wherein the constraints includes control invariant constraints that ensure satisfaction of specification constraints for the admissible values of the parameters of the model and admissible control inputs to the machine, and wherein the constraints include a terminal constraint specifying a region for a terminal state at the end of the prediction horizon, wherein the processor determines the region of the terminal constraint using a backward reachable method that iteratively determines, starting from a current set of feasible states of the machine, a backward reachable set of the states within the current set that after a one step evolution of dynamics of the machine is inside the current set until the backward reachable set of states of the machine is enclosed by the current set.

15. The system of claim 13, wherein the prediction horizon includes a sequence of prediction instants, such that the optimization produces control input for each prediction instant forming a sequence of control inputs for the prediction horizon, wherein estimates of the parameter of the model for determining the sequence of the control inputs changes between at least some prediction instants.

16. The system of claim 15, wherein the number of time instants is defined by a length of the prediction horizon.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1A is a block diagram of a controller for controlling an operation of a machine according to one embodiment of the invention;

(2) FIG. 1B is a block diagram of a general structure of the controller of FIG. 1A according to one embodiment of the invention;

(3) FIG. 1C is a block diagram of a method for controlling the operation of a machine according to one embodiment of the invention;

(4) FIG. 2A is a block diagram of the controller modules according to one embodiment of the invention;

(5) FIG. 2B is a block diagram of the module of the controller that constructs the prediction model according to one embodiment of the invention;

(6) FIG. 2C is a block diagram of the module of the controller that adjusts the cost function according to one embodiment of the invention;

(7) FIG. 2D is a block diagram of the module of the controller that constructs the constraints according to one embodiment of the invention;

(8) FIG. 2E is a block diagram of the module of the controller that computes the control input according to one embodiment of the invention;

(9) FIG. 2F is a flow chart of the method executed by the modules of the controller according to one embodiment of the invention;

(10) FIG. 2G is a flow chart of the method executed by the modules of the controller according to another embodiment of the invention;

(11) FIG. 3A is a schematic of an exemplar polytope arranged in the coordinate system visualizing the principles of some embodiments of the invention;

(12) FIG. 3B is a schematic of effect of the uncertainties of the parameters of the model on a state of the machine;

(13) FIG. 4A is a schematic of an example of a two-dimensional projection of the control-invariant subset of states defined by various constraints on the operation of the machine according to some embodiments of the invention;

(14) FIG. 4B is a schematic showing the relation between the feasible region of the states of the machine and control invariant subsets determined according to some embodiments of the invention;

(15) FIG. 5 is a block diagram of a backward-reachable region computation for determining the control invariant subset according to one embodiment of the invention;

(16) FIG. 6 is a block diagram of an exemplar implementation of determining the previous set of states according to one embodiment of the invention;

(17) FIG. 7 is a block diagram of a method for determining the region associated with the terminal constraint

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

(18) FIG. 1A shows a block diagram of a control system 101 for controlling an operation of a machine 102. The machine 102 is a device whose operation changes values of physical quantities such as positions, velocities, currents, and temperatures in response to commands. As used herein, the operation of the machine determines a motion of the machine that changes such quantities. The control system receives a desired motion 103 for the machine, such as a desired trajectory or target point for some of the quantities, and controls the machine via control inputs 104. The control inputs can include commands to change parameters of the operation of the machine or can include actual values of the parameters such as voltages, pressures, torques, forces that affect the machine motion resulting in the generation of output quantities 105 for the machine that are desired to match the desired motion 103. In some embodiments of the invention, the machine includes at least one parameter of uncertainty. For example, a model of an arm of a robot moving an object can include an uncertainty about a mass of the object carried by the arm. A model for the movement of a train can include an uncertainty about a friction of the wheels with the rails in current weather conditions.

(19) The control system 101 receives information 109 about the machine motion, from sensors, hardware, or software connected directly or remotely to the machine. The information 109 includes a state of the machine. The control system also receives the information about current estimates 108 of uncertain machine parameters from an estimation module 106. The control system 101 uses the state and the parameter estimates for the selection of the control inputs 104 which determine the motion of the machine and the output of the machine 105. The quantities 105, the control inputs 104 or a combination thereof, can be requested to remain in some pre-defined ranges according to constraints on the operation of the machine.

(20) FIG. 1B shows a general structure of the control system 101 according one embodiment of the invention. The control system 101 includes at least one processor 130 for executing modules of the controller. The processor 130 is operatively connected to a memory 120 for storing the model 112, the cost function 115 and the constraints 114 of the machine 101. It is an objective of some embodiments of the invention to determine the control inputs 104 using an adjustable cost function 115 and a model of the machine with uncertainty 112 subject to the constraints 114. In some embodiments, the model of the machine and the adjustable cost function are updated 116 based on current parameter estimates 108 during the operation of the machine. Alternatively, some embodiments update the model when the machine is not operating.

(21) Some embodiments of the invention are based on recognition that when the model 112 of the machine is updated, the cost function 115 also needs to be updated to achieve the convergence of the movement of the machine into a desired state. However such an update is generally computationally expensive. For example, such an update can require advanced numerical method running on powerful processors with large amounts of memory to find the solution of linear matrix inequalities or Riccati equations.

(22) FIG. 1C shows a block diagram of a method for controlling the operation of a machine according to one embodiment of the invention. The embodiment is based on realization that the cost function can be defined to include a parameter of the model having an uncertainty. In such a manner, when the parameter of the model is updated, such an update results in updating the cost function as well, which simplify the computational requirement for updating the cost function. The steps of the method can be performed by a processor such as the processor 130.

(23) The method iteratively controls 150 the operation of the machine with control inputs determined from an optimization of the cost function 115 along a prediction horizon subject to constraints on the control inputs. The optimization is performed using the model 112 according to the principles of the MPC.

(24) In this embodiment, the cost function 112 includes at least one parameter 155 of the model 112. The embodiment updates 160 the parameter 155 for at least some iterations of the controlling. In such a manner the update 160 result in updating both the model 112 and the cost function 115.

(25) FIG. 2A shows a block diagram of different modules of the control system 101 according to one embodiment of the invention. A model adjustment module 210 updates at least one parameter of the model 112 having the uncertainty based on current parameter estimates 108. Similarly, a cost function adjustment module 202 updates the cost function 115 to guarantee stability and convergence of the MPC. The outputs from a constraints and prediction horizon module 203 are used with the updated model 241 outputted by the module 201, the updated cost function 242 outputted by the module 202 and the machine state 109 to determine 204 the control input 104 to the machine

(26) FIG. 2B shows a block diagram of the model adjustment module 201. The current estimate of the parameters 108 is placed into a delay buffer 211 with length 218 equals to the length of the MPC prediction horizon and the oldest value 219 currently within the delay buffer is removed. The model of the machine with uncertainty is represented as a collection of models 212. The values of the current and past estimates of the parameters of the model stored in the buffer are used together with the collection of models to construct 213 the current MPC prediction model 241.

(27) According to some embodiments of the invention, the uncertainties of the parameters of the model are bounded. Those bounds on the uncertainties define ranges of admissible values of the parameters of the model, as explained in more details below. The admissible values form a collection of models, and some embodiments determine 213 the parameters of the model as a convex combination of the collection of the models.

(28) Thus, the collection of models is obtained as those models that represent all possible values of the uncertainties. For instance the collection of models can be obtained by considering the minimal and maximal value for the uncertain parameters and obtaining the models corresponding to all the admissible combinations of minimal and maximal parameters. A model for the machine with parameters within the range of the uncertain parameters can then be obtained.

(29) For example, each estimate of the parameter can be used together with the collection of models to obtain a model for the system valid at a specific time instant. Hence, by combining each value of the parameter estimates in the delay buffer with the collection of models, a sequence of models for the machine, where the sequence has length equal to that of the delay buffer, is obtained.

(30) In one embodiment of the invention, when new value of the parameter of the model is determined, the entire set of control inputs is determined using the model with updated parameter. This embodiment is simpler to implement and take advantage of the more accurate value of the parameter. However, alternative embodiments are based on recognition that such an update of the entire sequence of control inputs provides discontinuity with the sequence of control inputs previously determined with the old value of the parameter of the model.

(31) To that end, when the parameters of the model are updated, instead of predetermining the entire set of parameters in the buffer, some embodiments only remove the oldest 219 estimate of the parameters and add the new 218 estimate of the parameters. For instance considering an abstract representation of the buffer 211 in terms of sequential memory cells, when a new estimate is to be inserted, the parameter estimates currently in the buffer are copied to the memory cell to the right, when one exists. Thus, the oldest estimate placed on the rightmost cell 219 is overwritten, and hence disappears. The new estimate is placed on the leftmost memory cell 217, thus entering the buffer.

(32) This abstract representation of the buffer 211 in terms of sequential memory cells can correspond to different physical representation of memory cell in the memory 120 of the controller 101. Such an update leads to the situation when the buffer 211 stores different values of the parameters and thus the model 112 used for prediction is time varying along the steps of the prediction horizon, even if the uncertain parameter has a constant fixed value. However, this makes possible to update only a single component of the cost function by using only simple computations and still obtain stability.

(33) For example, the prediction horizon includes a sequence of prediction instants, such that the optimization produces control input for each prediction instant forming a sequence of control inputs for the prediction horizon. According to some embodiments, the model for determining the sequence of the control inputs changes between at least some prediction instants.

(34) For example, one embodiment determines, for each prediction instant, a prediction model using the parameter determined a number of time instants before the prediction instant to produce a sequence of prediction models for the prediction horizon and jointly determines the sequence of control inputs using the sequence of prediction models.

(35) In addition, some embodiments are based on realization that the proportionality between the estimation error and the control error and the convergence of the machine to the desired value are achieved if the machine parameters in the prediction model of the MPC are adjusted with a delay equal to the MPC prediction horizon, with respect to the time instant in which they are obtained from the estimation module. To that end, in one embodiment, the number of time instants is defined by a length of the prediction horizon.

(36) FIG. 2C shows a block diagram of a cost function adjustment module 202 according to one embodiment. The current cost function is computed 222 from the current parameter estimate 108 and a cost function components 221 representing the adjustable cost function 115. For example, in some implementations of this embodiment, the cost function includes a stabilization component 223 for directing a movement of the machine to a target state and a performance component 224 for optimizing the movement of the machine until the target state. In those implementations, the parameter 155 of the model includes a function of the stabilization component. For example, the stabilization component can include a sum of stabilization matrices weighted by an estimate of parameters of the model. For example, in addition to determining the stabilization matrices, one or more stabilizing gains can also be determined.

(37) FIG. 2D shows the details of the computation of the constraints and prediction horizon of MPC. From the specification constraints 231 the ensemble of cost function components, the preferred prediction horizon 232 and the ensemble of models a set of MPC constraints and actual prediction horizon for the MPC model are computed 233.

(38) FIG. 2E shows the details of the computation of the control input. From the adjusted model 241, adjusted cost function 242, prediction horizon 243, MPC constraints 244, and current machine state 245, a numerical problem is constructed 246, solved 247 by a numerical algorithm, and from the solution, the control input is selected 248.

(39) FIG. 2F shows a flowchart of the operation of the method in some embodiments of this invention. First the adjustable controller is designed 250 by determining collection of models 212, the cost function components 221, the constraints and the preferred horizon 232, and the parameters loaded in the control system memory 120, including the computed MPC constraints and prediction horizon 203. Then, the controlling of the machine begins. At each iteration, a new parameter estimate is determined 251, the current machine model is adjusted 252, the cost function is adjusted 253 and the machine state is determined 254. Then, the control input is computed 255, by first determining 256 a sequence of predicted inputs that satisfies the constraints and optimizes the cost function and then selecting 257 the initial part of such sequence. Then, the controller waits for next measurement 258.

(40) FIG. 2G describes a flowchart of the operation of the method in other embodiments of this invention. The operations are the same, however the adjustment of model and cost function occurs only once, before the selection of the first input, but after the controller memory has been programmed.

(41) Adjustment of Model

(42) Some embodiments of the invention are based on the realization that even if the true values of the parameters of the models are uncertain, the uncertainties of the parameters of the model are within a known range. For example, the mass of the train can be within a range of values for an empty train and fully loaded train. The ensemble of models 212 represents the models obtained for all the possible values of the uncertain parameters of the model.

(43) FIG. 3A shows an exemplar polytope 310 arranged in a coordinate system 320 visualizing the principles of the above realization. The coordinate system 320 is usually a low-dimensional system of parameters of the model having the uncertainties. For example, if the model includes four parameters with uncertainties, the coordinate system 320 is four dimensional coordinate system with an origin 325 and one dimension 322, 324, 326, and 328 for each of the parameters such that new values of the combination of parameters, e.g., values 333, 335, 337 are selected within the polytope 310.

(44) For example, a function that describes the motion of the machine at equally time-spaced sampling instants by a system of difference equations for a given machine state, input, and uncertainty d is
i. x(k+1)=(x(k),u(k),d(k))
y(k)=Cx(k)(1)
where k is the index of the sampling instant, is a generic function, C is a matrix, x is the machine state 106, u is the machine input 104 and y is the machine output 105. It is realized that in a number of situations the uncertainties are bounded within the range D, the states are bounded within the range , and the control inputs are bounded within the range U
ii. x,u,dD.(2)

(45) For instance the range of values D can be determined by the minimum and maximum value of the machine mass. Similarly, the ranges of the inputs and states can be determined by minimum and maximum values of velocities, positions, voltages, or torques.

(46) Some embodiments, based on and D, construct a set of matrices
{(A.sub.i,B.sub.i)}.sub.i=1.sup.l(3)
that results in the polytopic difference inclusion
i. co{A.sub.ix+B.sub.iu}.sub.i=1.sup.l(4)
such that for any given x and u,
ii. co{A.sub.ix+B.sub.iu}.sub.i=1.sup.l(x,u,d),dD,(5)
where co indicates the convex hull of the set.

(47) The polytopic difference inclusion (4) is described by the collection of models 212
i. {A.sub.ix+B.sub.iu}.sub.i=1.sup.l(4A)

(48) The linear models in Equation (4a) can be determined, for instance, by taking the maximum and minimum of the parameters that form vector d allowed by D, and/or of their combinations. Also, (4) can be constructed such that the desired value of the machine motion is placed at x=0.

(49) FIG. 3B shows a schematic of effect of the uncertainties of the parameters of the model determined by Equation (4) in estimating the state of the machine Given the current state x 301 of the machine and a control input, the next state of the machine can vary within a set 302 for different possible values of d of the parameters of the model of the machine. For instance, depending on specific values of the vector d, the state 302 can be transition to different states 310 within the set 302 for the same value of the control input.

(50) In some embodiments of the invention, the collection of models is designed such that the models for all values of the uncertain parameters within their ranges is included into a convex combination 303 of linear model. For example, the model of Equation (4) over-approximates the set 302 by set 303 which is a bounded polyhedron, or polytope, where the vertices 304 are determined from
i. .sub.i=A.sub.ix+B.sub.iu, i=1 . . . l(6)
wherein l, (A.sub.i, B.sub.i), i=1, . . . , l are the number and the matrices, respectively, defined in Equation (3).

(51) Based on Equation (4) an input u in the range that satisfies all
i. A.sub.ix+B.sub.iu, i=1(7)
makes the set 303 entirely in the range of admissible states, 320, and hence also the possible next state of the machine 302 according to (1) is within 320. Thus, enforcing constraints on Equation (4) guarantees enforcement of the constraints on Equation (1). The embodiments that use the model determined by Equation (4) are computationally advantageous, because the model determined by Equation (4) includes a set of linear models, rather than uncertain nonlinear models described in Equation (1).

(52) Thus in general the uncertain system model can be written as

(53) i.

(54) $\begin{array}{cc}x\left(t+1\right)=\underset{i=1}{\overset{l}{.Math.}}{\left[\right]}_i\left(A_{i}x+B_{i}u\right)& \left(8\right)\end{array}$
where is unknown but satisfies
ii.

(55) $\begin{array}{cc}{l}^{},0{\left[\right]}_{i}1,\underset{i}{.Math.}{\left[\right]}_i^l=1,& \left(9\right)\end{array}$
such that the uncertain system model is the convex combination of the ensemble of models (7), and [a].sub.i denotes the i.sup.th component of vector a.

(56) The model for a specific realization of the uncertain parameters is written

(57) i.

(58) $\begin{array}{cc}\hat{A}x+\hat{B}u=\underset{i=1}{\overset{l}{.Math.}}{\left[\right]}_i\left(A_{i}x+B_{i}u\right),& \left(10\right)\end{array}$
thus, the estimation module 106 provides a current estimate of the parameters at time t by determining (t) and the corresponding estimated model (10).

(59) For instance, the module 106 can find (t) as the value that provide the best fit while considering a recent set of input and output or state data, possibly with some move reduction term

(60) $\left(t\right)=\arg \underset{}{\min }\underset{k=t-N_e}{\overset{t}{.Math.}}||x\left(k\right)-m\left(,k\right)||s.t.m\left(,k\right)={\left[\right]}_i{A}_{i}x\left(k-1\right)+\mathrm{Bu}\left(k-1\right)$

(61) However, the embodiments of this invention are not restricted to a specific estimation module, as they work for any estimation module as long as this provides a valid estimate (t) that satisfies (9).

(62) Adjustment of Cost Function

(63) The cost function of MPC controller encodes the desired priority of regulating the states of the machine to achieve the desired machine behavior while limiting the use of the machine actuation system so to save energy, minimize vibration, reduce wear of parts.

(64) Thus, for models such as (10) where is fixed and known when the controller is designed the MPC cost function can be written as

(65) $\begin{array}{cc}\underset{k=0}{\overset{N-1}{.Math.}}{x\left(k\right)}^{}\mathrm{Qx}\left(k\right)+{u\left(k\right)}^{}\mathrm{Ru}\left(k\right),& \left(11\right)\end{array}$
where Q, R are weighting matrices that encode the desired priority of regulating the states of the machine

(66) However, the cost function (11) that encodes the performance objective of the controlling of the machine needs to be modified to obtain guarantees of converging to the desired state and stability of the machine operation. For example, it is possible to add a stabilizing term to (11) which for models such as (8) results in
x(N)Px(N)(12)

(67) Thus, the complete cost function is the combination of the performance term (11) and the stability term (12)

(68) $\begin{array}{cc}{x\left(N\right)}^{}\mathrm{Px}\left(N\right)+\underset{k=0}{\overset{N-1}{.Math.}}{x\left(k\right)}^{}\mathrm{Qx}\left(k\right)+{u\left(k\right)}^{}\mathrm{Ru}\left(k\right)& \left(13\right)\end{array}$
where P is an appropriately designed weighting matrix for the terminal state at the end of the MPC prediction horizon.

(69) Some embodiment of the invention are based on realization that while the model changes due to changing estimates of the system parameters, the performance term (11) can remain unchanged, but the stabilizing component (13) needs to be adjusted to achieve stability of the updated model. Accordingly, some embodiments modify the stabilizing component to include a stabilizing term that depends on the parameter estimate that is
x(N)P((N))x(N)(14)
where P() is a parameter dependent weight matrix.

(70) According to some embodiments, the term P() for the collection of models (7) and the uncertain model (8), (9) includes convex combination of matrix weights obtained from the current estimate of the parameter ,

(71) $\begin{array}{cc}P\left(\right)=\underset{i=0}{\overset{l}{.Math.}}{\left[\right]}_i{P}_i& \left(15\right)\end{array}$

(72) It is further realized that for the case where B.sub.i=B, ii=1, . . . , l, the matrix weights in (15) can be designed by solving the linear matrix inequality

(73) $\begin{array}{cc}\left[\begin{array}{cccc}{G}_i+G_i^{}-S_i& {\left(A_i{G}_i+{\mathrm{BE}}_i\right)}^{}& E_i^{}& G_i^{}\\ \left(A_i{G}_i+{\mathrm{BE}}_i\right)& S_j& 0& 0\\ E_i& 0& R^{-1}& 0\\ G_i& 0& 0& Q^{-1}\end{array}\right]>0,i,j{}_{\left[1,l\right]}.& \left(16\right)\end{array}$
for variables G.sub.i, S.sub.i.sup.nn, S.sub.i>0, E.sub.i.sup.mn, i.sub.[1, l], and then computing
P.sub.i=S.sub.i.sup.1, K.sub.i=E.sub.iG.sub.i.sup.1, i.sub.[1,l],(16b)
to construct the matrices in (15).

(74) Instead, the liner matrix inequality can be represented as

(75) $\begin{array}{cc}\left[\begin{array}{cccc}G+G^{}-S_i& {\left(A_{i}G+\mathrm{BE}\right)}^{}& E^{}& G^{}\\ \left(A_{i}G+\mathrm{BE}\right)& S_j& 0& 0\\ E& 0& R^{-1}& 0\\ G& 0& 0& Q^{-1}\end{array}\right]>0,i,j{}_{\left[1,l\right]}.& \left(17\right)\end{array}$

(76) S.sub.i.sup.nn, S.sub.i>0, E, G.sup.mn, is solved and
P.sub.i=S.sub.i.sup.1, i.sub.[1,l], K=EG.sup.1.(17b)

(77) The complexity of linear matrix inequalities can grow for large-scale systems. However, some embodiment's exploit symmetries of the model of the machine to reduce computational cost. A symmetry is a pair of invertible matrices and such that for all i.sub.[1, l] there exists j.sub.[1, l] that satisfy
A.sub.i.sup.1=A.sub.j
B.sup.1=B
X=X
U=U.

(78) Some embodiments of the invention are based on the realization that because the symmetry relates extreme systems models A.sub.i, the extreme cost functions and controllers can be related. Accordingly, one embodiment determines at least part of the stabilizing term (14) of the cost function (13) using a symmetry between values of the parameter of the model to reduce the computations in solving (16), (17).

(79) For example, some embodiments solve the equations (16), (17) only for the models in the collection that are not related by symmetry, thus reducing the number of variables in (16), (17) and thus making them easier to solve. For example, denote by .sup.1.Math..sub.[1, l] a subset of the models A.sub.i in the collection of models that are not related by symmetry. Also, denote by .sup.2.Math..sup.2.sub.[1, l] a subset of pairs (i, j) of models in the collection that are not related by symmetry. Some embodiments determine variables S.sub.i, G.sub.i, and E.sub.i only once for each symmetrically unrelated model i.sup.1 in the collection. Also, some embodiments solve the equations (15) only once for each symmetrically unrelated pair (i,j).sup.2.

(80) For example, one embodiment solves

(81) 0$\left[\begin{array}{cccc}{G}_i+G_i^{}-S_i& {\left(A_i{G}_i+{\mathrm{BE}}_i\right)}^{}& E_i^{}& G_i^{}\\ \left(A_i{G}_i+{\mathrm{BE}}_i\right)& S_k{}^{-1}& 0& 0\\ E_i& 0& R^{-1}& 0\\ G_i& 0& 0& Q^{-1}\end{array}\right]>0$
for each (i, j).sup.2 where i.sup.1 and S.sub.k.sup.1=S.sub.j. Then the embodiment determines
E.sub.k=E.sub.i.sup.1
G.sub.k=G.sub.i.sup.1
S.sub.k=S.sub.i.sup.1
for each k.Math..sup.1 where A.sub.k=A.sub.i.sup.1 for i.sup.1.

(82) Constraints Computation

(83) Some embodiments of the invention are based on realization that in order to guarantee that the constraints on the machine motion are satisfied during the operation of the machine for all variation of the values of the parameters of the model within the predetermined range, a subset of the feasible region of the states can be determined, such that whenever the state of the machine is in such a subset, there exists at least one admissible control input such that the state remains in the subset for all possible values of the parameters and the disturbances. By remaining in the subset, all the machine motion constraints are satisfied because the subset is included in the feasible region .

(84) Some embodiments of the invention determine and use such constraints on the control inputs, referred herein as control-invariant constraints on the control inputs, to determine a control input that maintains the state of the machine in such a subset of the feasible region of the states, referred herein as a control-invariant subset of states. By the construction, a control input satisfying the control-invariant constraints maintains a state of the machine in a control-invariant subset of states satisfying constraints on the operation of the machine for all values of the parameter of the model within the range defined by the uncertainty model.

(85) FIG. 4A shows an example of a two-dimensional projection of the control-invariant subset of states 410 defined by various constraints on the operation of the machine according to embodiments of the invention. Typically, the feasible region is a multi-dimensional polytope determined by hyperplanes, which are represented by linear inequalities, along multiple dimensions corresponding to the constraints on the operation.

(86) Due to the nature of receding horizon control, the existence of a solution for a certain horizon does not by itself guarantees the existence of the solution for a subsequent horizon. For example, the state of the machine and a state of the reference trajectory 420 can be optimal and feasible for one iteration, but all control actions 421-424 that controller is allowed to take during the next iteration can bring a state of the machine outside of the feasible region 410.

(87) Some embodiments of the invention are based on yet another realization that it is possible to select a control-invariant subset 415 of the feasible region, such that from any state of the machine within that control-invariant subset, there is a control input maintaining the state of the machine within the subset for the known future states of the reference trajectory or for all admissible future states of the reference trajectory. For example, for any state such as a state 430 within the subset 415 and within all possible control inputs 431-434 that the controller can execute, there is at least one control input 434 that maintains the state of the machine and reference within the subset 415. In this case, the subset 415 is a control invariant subset.

(88) Some embodiments determine the control invariant subset of state for the model of Equation (4), such that for each state within the control invariant subset there is at least one control action maintaining the state of the machine within the control invariant subset for all possible states values according to Equation (5).

(89) FIG. 4B shows the relation between the feasible region 401 of the states of the machine, the largest control invariant subset 402 of the state for the model with uncertainties according to Equation (1) within . In some embodiments the control invariant subset 402 is reduced to form control invariant subset 403, which is less than the subset 402, but has a simpler shape formed by linear equations. Because the subset 402 is determined by nonlinear equations, the shape of the subset 402 can be non-convex. Thus, computation of 403 is simpler, more efficient and advantageous for control purposes.

(90) FIG. 5 shows a block diagram of a backward-reachable region computation for determining the control invariant subset 403 starting from the feasible region , 401 according to one embodiment of the invention. The backward-reachable region computation determines the robust control invariant subset .sub.x and also the control invariant set .sub.u(x), which determines for any x within .sub.x the set of inputs u in that can be applied so that all possible next states according to (4) are inside .sub.x.

(91) The backward-reachable region computation initializes 501 a current set .sub.c to the feasible set and determines 502 a previous set of states .sub.p as a subset of the current set .sub.c such that for all states x in .sub.p there exists an input u in such that for all the possible values of the parameters p in P, the updated state is in the current set .sub.c.

(92) If 503 the previous set .sub.p is empty 504, correct operation of the controller cannot be guaranteed, which means that the set P of possible values of the parameters should be reduced in size, possibly by changing the design or objective of the operation of the machine. If 505 the current set and the previous set are equal, that is also 506 the set .sub.x otherwise, the previous set is assigned 507 to be the current set and the computation iterates 508 again. When the set .sub.x is found the last computed set of state-input couples is the robust admissible input set .sub.u(x) for all x within .sub.x.

(93) FIG. 6 shows a block diagram of an exemplar implementation of determining 502 the previous set of states according to one embodiment. The embodiment identifies 601 the state-input couple that generates an updated state that is in the current set for all the values of the parameters, and projects 602 the state input couples into state values, i.e., the embodiment identifies the states that belongs to at least one of such state-input couples.

(94) A method of computation of the couples of states-inputs 601 is based on determining the couples (x, u) such that the updated state is inside the current set for all the models in the ensemble, i.e., A.sub.ix+B.sub.iu.sub.c, i=1, . . . , l.

(95) In addition to the specification constraints on the operation of the machine, some embodiments add additional constraints to the MPC in order for the machine state to converge to the desired value. For example, such an additional constraint is a terminal constraint, such as
x(N).sub.N,(18)
which determines a region for the terminal state at the end of the MPC prediction horizon.

(96) In some embodiments, such terminal constraint is determined according to the specification constraints and the model of the machine. However, such a determination is difficult when the model of the machine is uncertain. Thus, an alternative embodiment of the invention determines the terminal constraint using the specification constraints 231, the collection of models 211 and the cost function components 221.

(97) For example, one embodiment determines the terminal constraints as a set
.sub.N={x:(A.sub.i+B.sub.iK.sub.i).sup.tx,K.sub.ixU, i,j.sub.[1,l],t[0,)},(19)
which is the set of initial states where the machine controlled from the control gains K.sub.i, i=1, . . . , l or K in (16b) or (17b) obtained by solving (16) or (17), respectively, satisfies the constraints.

(98) Some embodiments determine the set of the terminal constraint as the fixed point of a sequence of backward reachable sets by
.sup.(0)={x:(x,K.sub.ix).sub.xu,i.sub.[1,l]}
.sup.(h+1)={x:(A.sub.i+BK.sub.i)x.sup.(h),i.sub.[1,l]}.sup.(h),(20)
where .sub.xu.Math.U is any subset of feasible states and inputs.

(99) FIG. 7 shows a block diagram of a backward reachable method for determining a region corresponding to the terminal constraint according to one embodiments of the invention. The backward reachable method iteratively determines, starting from a current set of feasible states of the machine, a backward reachable set of the states within the current set that after a one step evolution of dynamics of the machine is inside the current set until the backward reachable set of states of the machine is enclosed by the current set.

(100) For example, a current set is initialized 701 as a set of feasible states where the input generated by control gains K.sub.i, i=1, . . . , l or K in (16b) or (17b) obtained by solving (16) or (17), respectively, satisfy input constraints. From the current set, a backward reachable set is computed 702 to determine the states within the current set that after a one step evolution of the system dynamics are inside the current set. If 703 the backward reachable set contains or equals the current set, then 704 the current set is the maximal control invariant set. Otherwise 705, the backward reachable set intersects with the current set, and used as the new current set, in a repetition 706 of the method.

(101) In some embodiments, due to the presence of invariant constraints, the backward reachable method in FIG. 7 is initialized according to
.sub.xu={(x,u),u.sub.u(x)}.Math.U.

(102) However, in order to guarantee that the terminal constraint (18) can always be satisfied, some embodiments modify the prediction horizon. Such an adjustment helps to ensure that the terminal constraints can be satisfied, that is, that the region corresponding to the terminal set can be reached by any feasible initial state within the length of the prediction horizon.

(103) To that end, some embodiments determine a minimal length of the prediction horizon for which there is a solution to the optimization of the cost function for all states satisfying the control invariant constraints. Next, the embodiments determine a length of the prediction horizon as largest between a preferred length and the minimal length.

(104) For example, some embodiments determine the minimum prediction horizon N.sub.m such that the equation (3) can always be satisfied for x(0).sub.x and if the preferred horizon N.sub.p 232 is less than N.sub.m, the embodiments set N=N.sub.m, otherwise the embodiments set N=N.sub.p. The computation of the minimum horizon can be obtained by determining the minimum index N.sub.m such that the set .sup.(N.sup.m.sup.) determined from the sequence of intersected backward reachable sets according to

(105) $\begin{array}{cc}\begin{array}{c}{\mathrm{��}}^{\left(0\right)}={\mathrm{��}}_N,\\ {\mathrm{��}}_i^{\left(h+1\right)}=\left\{x\mathrm{��}\text{:}u\mathrm{��},A_{i}x+\mathrm{Bu}{\mathrm{��}}^{\left(h\right)}\right\},\\ {\mathrm{��}}^{\left(h+1\right)}=\underset{i=1}{\overset{l}{.Math.}}{\mathrm{��}}_i^{\left(h+1\right)}.\end{array}& \left(21\right)\end{array}$
is such that .sup.(N.sup.m.sup.)=.sub.x.

(106) Some embodiments rake advantage from the realization that if the terminal set .sub.N is control invariant, that is, it is invariant for any combination of models in the collection and gains K.sub.j obtained from equations (16), (17), then a control invariant .sub.x that also guarantees that the constraint (18) is always satisfied can be computed for the preferred prediction horizon.

(107) In one embodiment, this computation is obtained by constructing N.sub.p-steps backward reachable set by performing exactly N.sub.p iterations of the equations
.sup.(0)=.sub.N,
.sup.(h+1)={x: uU,A.sub.ix+Bu.sup.(h),i.sub.[1,l]}.(33b)

(108) In this case the MPC horizon can be set equal to the preferred horizon set N=N.sub.p and .sub.x=.sup.(N).

(109) For instance if B.sub.i=B, ii=1, . . . , l and (17) is solved, then

(110) .sub.N computed from (20) is robust invariant and the method in (33b) can be used.

(111) Prediction Model

(112) Some embodiments of the invention are based on realization that since the MPC controller needs to predict the evolution of the model along its prediction horizon N, the estimate of the model at time t cannot be used for all prediction steps. Instead, the current estimate at time t, (t) is placed into the delay buffer 211, that has the length of the MPC prediction horizon N and that includes the most recent N values of the estimates . Then, the construction of the numerical problem 246 builds the MPC prediction model from all the past estimates placed on the delay buffer.

(113) Thus, the prediction model k steps in the future with respect to control time t is obtained as

(114) $\begin{array}{cc}{x}_{k+1|t}=\underset{i=1}{\overset{l}{.Math.}}{\left[\left(t-N+k\right)\right]}_i\left(A_i{x}_{k|t}+B_i{u}_{k|t}\right)& \left(22\right)\end{array}$
for all k[0, N].

(115) Thus, in some embodiments the model used for determining the control inputs for the prediction horizon varies for different steps along the prediction horizon. For example, the model of equation (22) is a linear time varying systems constructed from the collection of the models and the N most recent estimates of the parameter, where at each step along the prediction horizon, the parameter estimate obtained N steps before is used.

(116) Accordingly, the stabilization term (14) of the cost function is adjusted from (t), by setting (N)=(t) in (14) and evaluating (14) by (15).

(117) For the case where the adjustment of the parameters occurs only before of the operation of the machine, in this invention it is realized that it is enough to set
(tN+k)={circumflex over ()},t,N(23)
where {circumflex over ()} is the estimate obtained at the moment the estimation is terminated, and use such value for both model adjustment (22) and cost function adjustment (14).

(118) Determining Control Inputs

(119) In some embodiments, from the prediction model, cost function, MPC constraints, and the current state x(t) the MPC controller constructs a finite horizon optimal control problem

(120) $\begin{array}{cc}{\mathrm{��}}_{{}_t^N}^{\mathrm{MPC}}\left(x\left(t\right)\right)=\underset{U_t}{\min }{x}_{N|t}^{}P\left({}_{N|t}\right)x_{N|t}+\underset{k=0}{\overset{N-1}{.Math.}}{x}_{k|t}^{}{\mathrm{Qx}}_{k|t}+u_{k|t}^{}{\mathrm{Ru}}_{k|t}s.t.x_{k+1|t}=\underset{i=1}{\overset{l}{.Math.}}{\left[\left(t-N+k\right)\right]}_i{A}_i{x}_{k|t}+{\mathrm{Bu}}_{k|t}u\mathrm{��},x\mathrm{��}{u}_{k|t}{\mathrm{��}}_u\left(x_{k|t}\right)x_{N|t}{\mathrm{��}}_N{x}_{0|t}=x\left(t\right),& \left(24\right)\end{array}$
where U.sub.t=[u.sub.0|t . . . u.sub.N1|t] is the sequence of predicted inputs, .sub.t.sup.N=[(tN) . . . (t)], is the sequence of N most recent predictions on system parameters and .sub..sub.N.sup.MPC(x) is the nonnegative value function that holds 0 value only if the state of the machine is at the desired value.

(121) The problem is further converted into a numerical problem, which is a quadratic program as

(122) $\begin{array}{cc}\underset{U_t}{\min }\frac{1}{2}{U}_t^{}{Q}_p{U}_t+C_p^{}{U}_{t}s.t.G_p{U}_t{W}_p.& \left(25\right)\end{array}$

(123) The solution of the quadratic problem (25) provides the optimal input sequence for the current state and parameter estimate U*.sub.t=[u*.sub.0|t . . . u*.sub.N1|t]. The control input for the machine at current time t is then selected as the first element of the optimal control sequence, e.g.,
u(t)=u*.sub.0|t.(26)

(124) The capability of computing the solution to (24) is based on the fact that the obtained numerical problem (26) is convex whenever the cost function matrices Q, R are positive definite, for all value of machine state and all admissible estimates of the machine parameters. Also, the satisfaction of the constraints in presence of parameter estimation error is due to the construction of .sub.x and .sub.N. The convergence of the machine to the desired state is due to the construction of the adjusted cost function (13) and terminal constraint (18) according to the embodiments, which guarantees that
.sub..sub.t+1.sub.N.sup.MPC(x(t+1)).sub..sub.t.sub.N.sup.MPC(x(t)).sub.min(Q)x(t).sup.2+.sub.ISS{circumflex over ()}(tN),(27)
where .sub.min(Q) is the minimum eigenvalue of Q, which is greater than 0 since Q is positive definite, and ISS is the sensitivity gain from the estimation error.

(125) The above-described embodiments of the present invention can be implemented in any of numerous ways. For example, the embodiments may be implemented using hardware, software or a combination thereof. When implemented in software, the software code can be executed on any suitable processor or collection of processors, whether provided in a single computer or distributed among multiple computers. Such processors may be implemented as integrated circuits, with one or more processors in an integrated circuit component. Though, a processor may be implemented using circuitry in any suitable format.

(126) Further, it should be appreciated that a computer may be embodied in any of a number of forms, such as a rack-mounted computer, a desktop computer, a laptop computer, minicomputer, or a tablet computer. Such computers may be interconnected by one or more networks in any suitable form, including as a local area network or a wide area network, such as an enterprise network or the Internet. Such networks may be based on any suitable technology and may operate according to any suitable protocol and may include wireless networks, wired networks or fiber optic networks.

(127) Also, the various methods or processes outlined herein may be coded as software that is executable on one or more processors that employ any one of a variety of operating systems or platforms. Additionally, such software may be written using any of a number of suitable programming languages and/or programming or scripting tools.

(128) Also, the embodiments of the invention may be embodied as a method, of which an example has been provided. The steps performed as part of the method may be ordered in any suitable way. Accordingly, embodiments may be constructed in which acts are performed in an order different than illustrated, which may include performing some acts simultaneously, even though shown as sequential acts in illustrative embodiments.

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