DETAILED DESCRIPTION OF THE INVENTION

[0030] Referring to FIG. 2 , there is shown in block diagram form a implementation of Model Predictive Control (MPC) block 50 from the control hierarchy of FIG. 1 . As noted above and as is apparent from FIG. 2 , MPC block 50 is divided into a steady-state calculation and a dynamic calculation. A “plant” is also represented by block 60 in FIG. 2 . As noted above, as used herein the term “plant” is intended to refer to any of a variety of systems, such as chemical processing facilities, oil refining facilities, and the like, to which those of ordinary skill in the art would appreciate that model predictive control technology may be advantageously applied. The exact nature of plant 60 is not important for the purposes of the present disclosure.

[0031] In the interest of clarity, not all features of actual implementations of an MPC controller are described in this specification. It will of course be appreciated that in the development of any such actual implementation, as in any such project, numerous engineering and programming decisions must be made to achieve the developers' specific goals and subgoals (e.g., compliance with system- and business-related constraints), which will vary from one implementation to another. Moreover, attention will necessarily be paid to proper engineering and programming practices for the environment in question. It will be appreciated that such a development effort might be complex and time-consuming, but would nevertheless be a routine undertaking for those of ordinary skill in the field of controller system design having the benefit of this disclosure.

[0032] The steady-state target calculation for linear MPC, represented by block 62 in FIG. 2 , takes the form of a linear program (“LP”). The majority of academic literature focuses on the dynamic problem, represented by block 64 in FIG. 2 , and not the steady-state problem; it is left to industrial practitioners to investigate issues associated with the steady-state LP. Because most models used in industry are linear, with linear economics driving the controller, the result is an LP. (As noted above, it is contemplated that the present invention may be advantageously practiced in application involving non-linear (e.g., quadratic or higher-order) objectives.) Briefly, the overall system depicted in FIG. 2 , comprising model predictive controller 50 and plant 60 , operates as follows: MPC 50 performs dynamic and steady-state calculations to generate control signals reflecting optimal “inputs” u* to plant 60 . As represented in FIG. 2 , inputs u* are carried to plant 60 on line 66 . Thus, line 66 in FIG. 2 is representative of what would be, in a “real-world” application of MPC technology, a plurality of electrical control signals applied to controlled components (valves, for example) to control a plurality of controlled variables (pressures, flow rates, temperatures, for example) within the plant.

[0033] On the other hand, the plant 60 's operation is symbolically represented by a plurality of “outputs” y which are represented in FIG. 2 as being carried on a line 68 . thus, line 68 is representative of a plurality of electrical signals reflecting the operational status of plant 60 .

[0034] Note that in the system of FIG. 2 , output(s) y are fed back to be provided as inputs to steady-state target calculation block 62 . Steady-state target calculation block 62 operates to generate so-called “target” inputs and outputs u _s and y _s , respectively, as a function of plant output(s) y and as a function of an “objective” which is symbolically represented as being applied to calculation block 62 on line 72 . The target inputs u _s and y _s , are represented in FIG. 2 as being carried on a line 70 . Target inputs u _s are those inputs to plant 60 which, based on the calculations made by calculation block 62 , are expected based on the MPC modeling of plant 60 , to result in plant 60 operating to produce the target outputs y _s .

[0035] Those of ordinary skill in the art will appreciate that there may be certain real-world constraints or considerations which would make it undesirable for target inputs u _s to be applied directly to plant 60 . It is the function of dynamic regulator 64 to derive, as a function of target inputs and outputs u _s and y _s , “optimal” inputs u* to plant 60 .

[0036] A simple example is illustrative of what FIG. 2 represents. In this example, “plant” 60 is a water faucet. “Output” y is a value reflecting the rate of flow of water out of the faucet. The “objective” in this example may be to keep the output (i.e., flow rate) at a particular non-zero value. Given this objective, steady state target calculation block 62 derives target inputs and outputs u _s and y _s that, based on the MPC model, can be expected to achieve outputs y that accomplish the objective. The single “input” in this example would be a value reflecting the extent to which the faucet is turned on. Dynamic regulator 64 , in turn, calculates an optimal input u* as a function of target inputs and outputs u _s and y _s that, based on the dynamic calculations of regulator 64 , will result in target outputs y _s .

[0037] The water faucet example illustrates the different roles of steady-state control and dynamic control. Suppose, for example, that the objective for the hypothetical system is to keep the water flowing out of the faucet at a given rate of x cm ³ /second. Assuming that the system begins operation with the faucet turned off (i.e., y=0 at t=0), steady-state target calculation block 62 would generate a target input us that, based on its MPC model of the faucet, would result in a flow rate of x cm ³ /second. If this target input were applied directly to the faucet, this would require the faucet to instantly be opened to the extent necessary to achieve the objective flow rate. Those of ordinary skill in the art will appreciate that this may be neither desirable nor possible. Instead, dynamic regulator 64 operates to derive an “optimal” input u* based on the target inputs and outputs u _s and y _s provided on line 70 . In the example, the “optimal” input may reflect a more gradual turning on of the faucet. Over time, there may be expected to be a convergence of the “optimal” input u* and the target input u _s .

[0038] In the foregoing simplistic example, there was but a single input to plant 60 and a single output from plant 60 . It is to be understood, however, that in typical embodiments, there are likely to be a plurality of inputs and a plurality of outputs from the plant. In those cases, the values u and y are vectors.

[0039] The MPC system depicted in FIG. 2 represents a so-called “input/output” model of control technology. In the input/output model, the current output of the system (i.e., plant 50 ) is expressed as a function of past inputs to the system. Expressed mathematically, y _k =ƒ(u _k , u _k−1 , u _k−2 , . . . ). An alternative model of control technology is known as the “state space” model. In a state space model, the current state of the system is expressed as a function of the previous “state” of the system and the previous input to the system. The “state” of the system, represented by the variable x, is one or more values reflecting the current operational status of the system. The “state” of the system is not strictly speaking the same as the “output” of the system, although the two terms are not altogether dissimilar. There may be state variables which are not outputs from the system or are not even perceptible from outside the system. Expressed mathematically, in the state space model, the current state x _k+1 is some function F of the previous state x _k and the previous input u _k (that is, x _k+1 =F(x _k , u _k )), while the current output y _k is some function G of the current state x _k (that is, y _k =G(x _k )).

[0040] Turning now to FIG. 3 , there is shown an MPC system based on the state space model. The system represented in FIG. 3 includes an MPC controller 50 ′ which differs in some respects from controller 50 in FIG. 2 , although in the embodiment of FIG. 3 , those elements which are substantially the same as in that of FIG. 2 have retained identical reference numerals. The principal difference between controllers 50 and 50 ′ is the inclusion in controller 50 ′ of an estimator block 74 which receives as inputs both the output(s) y of plant 60 on line 68 , as well as the “optimal” input(s) u* to plant 60 on line 66 . Based on these inputs, estimator 74 calculates and/or estimates the current state of plant 60 and provides this information and the output on line 76 to steady-state target calculator 62 .

[0041] In the embodiment of FIG. 3 , steady-state target calculator 62 calculates target states x _s and/or target outputs y _s and target inputs u _s . As in the embodiment of FIG. 2 , dynamic regulator 64 uses the output of steady-state target calculator 62 to compute optimal inputs u* to plant 60 .

[0042] Notably, as will hereinafter become apparent to those of ordinary skill in the art, the present invention is independent of the model used, and may be advantageously practiced in the context of either the input/output model or the state space model. In particular, the present invention relates to the operation of the steady state calculation block 62 in the systems of FIGS. 2 and 3 .

[0043] Those of ordinary skill in the art will appreciate that any linear model, whether it be state-space, step response, impulse response, or other, can be cast—at steady-state—in the following form:

Δy=GΔu (1)

[0044] Here Δy=[Δy ₁ . . . Δy _n ] ^T ∈ ^{n _y} represents the change between the current steady-state output and the last, and Δu=[Δu ₁ . . . Δu _n ]∈ ^{n _u} represents the change between the current steady-state input and the last:

Δu=u−u _k−1 (2)

Δy=y−y _k−1 (3)

[0045] Here y=└y ₁ . . . y _{n y} ┘ ^T and u=└u ₁ . . . u _{n u} ┘ ^T are vectors containing the future steady-state outputs and inputs, respectively, and y _k−1 and u _k−1 are the previous measured values of y and u. G∈R ^{n _y xn _u} is the steady-state gain matrix for the system. The gain matrix G expresses the sensitivity of or relationship between changes in output to changes in input. The LP tries to find the optimal steady-state targets by minimizing or maximizing (i.e., extremizing) an objective function J _s which may represent, for example, an economic consideration associated with operation of the controlled system: 1 $\begin{array}{cc}{J}_s=\sum _{i=1}^{n_u}c_iu_i+\sum _{j=1}^{n_y}d_jy_j& \left(4\right)\end{array}$

[0046] while maintaining the relationship of Equation (1) above between the steady-state inputs and outputs, and respecting input and output constraints arising from process specifications:

+E,uns u _i ≦u _i ≦{overscore (u _i )} i=1, . . . , n _u (5)

+E,uns y _j ≦y _j ≦{overscore (y _j )} j=1, . . . , n _y (6)

[0047] and ensuring the resulting solution does not force the dynamic MPC calculation to become infeasible:

NΔ+E,uns u _i ≦Δu _i ≦NΔ{overscore (u _i )} i=1, . . . , n _u (7)

[0048] (As used herein, the term “feasible” as applied to system-manipulated (input) variables and system-controlled (output) variables means that the variables do not violate upper and lower constraints arising from predetermined process specifications, such as those set forth in equations (5) and (6) above. Conversely, “infeasible” refers to a violation of these constraints. Such constraints constitute bounds that limit the variables to some previously determined process operating region.)

[0049] In Equation (7), N refers to the control horizon of the dynamic MPC calculation, and Δ+E,uns u _i and Δ{overscore (u)} _i are the minimum and maximum bounds for the rate of change constraints in the dynamic MPC calculation. In Equations (5) and (6), +E,uns u _i and {overscore (u)} _i are minimum and maximum bounds for u _i , respectively, with similar notation employed for the outputs y _i . Equation (7) ensures the steady-state target is compatible with the rate of change or velocity constraints in the dynamic MPC calculation.

[0050] To avoid real-time infeasibilities in the LP, it is common to recast the hard output constraints +E,uns y _j ≦y _j ≦{overscore (y)} _j in Equation (7) as a soft constraint by adding a so-called “slack” variable ε=└ε ₁ . . . ε _{n y} ┘ ^T ≧0∈ ^{n _y} that allows for some amount of violation in the constraint:

−ε _j ++E,uns y _j ≦y _j ≦{overscore (y)} _j +ε _j j=1, . . . , n _y. (8)

[0051] The size of the violation is minimized by appending it to the objective function: 2 $\begin{array}{cc}{J}_S=\sum _i^{n_u}c_iu_i+\sum _j^{n_y}d_jy_j+\sum _j^{n_y}{\varepsilon }_j.& \left(9\right)\end{array}$

[0052] where c ₁ represents costs or weights of input variables u _i as they relate to the objective function, and d _i represents such costs or weights of output variables y _i .

[0053] Additionally, in industrial implementations, a bias may be introduced to the model as a result of output feedback. It is assumed that the difference between the model prediction and the measured output at the current time is due to a constant step disturbance at the output. While there are other possible types of disturbance models that can be incorporated into the MPC framework; the present disclosure will concentrate on output disturbance models. The model bias

b=└b ₁ . . . b _{n y} ┘ ^T ∈ ^{n _y} (10)

[0054] is based on a comparison between the current predicted output y and the current measured output ŷ=└ŷ ₁ . . . ŷ _{n y} ┘ ^T ∈ ^{n _y} :

b=ŷ−y. (11)

[0055] The bias is added to the model from Equation (1) above:

Δ y=GΔu+b. (12)

[0056] The problem is then: 3 $\begin{array}{cc}\underset{u,y,\epsilon }{\min }\sum _i^{n_u}c_iu_i+\sum _j^{n_y}d_jy_j+\sum _j^{n_y}{\varepsilon }_j& \left(13\right)\end{array}$

[0057] subject to 4 $\begin{array}{cc}\Delta y=G\Delta u+b\underset{\_}{u_i}\le u_i\le \stackrel{\_}{u_i}N\Delta \underset{\_}{u_i}\le \Delta u_i\le N\Delta \stackrel{\_}{u_i}-{\varepsilon }_J+\underset{\_}{y_J}\le y_j\le \stackrel{\_}{y_J}+{\varepsilon }_J0\le {\varepsilon }_j& \left(14\right)\end{array}$

[0058] with i=1, . . . , n _u and j=1, . . . , n _y . The LP is generally considered simpler than the dynamic MPC calculation (which takes the form of a quadratic program) and is commonly solved using the mainstay of linear programming, the simplex method. It is common to calculate the moves Δu and Δy, instead of the inputs and outputs directly. The final nominal steady-state target calculation can be expressed as: 5 $\begin{array}{cc}\underset{\Delta u,\Delta y,\varepsilon }{\min }c^T\Delta u+d^T\Delta y+e^T\varepsilon & \left(15\right)\end{array}$

[0059] subject to 6 $\begin{array}{cc}\Delta y=G\Delta u+bA_u\Delta u\le b_uA_y\Delta y\le b_y+\varepsilon \varepsilon \ge 0& \left(16\right)\end{array}$

[0060] where

c=└c ₁ . . . c _{n u} ┘ ^T ∈ ^{n _u} (17)

d=└d ₁ . . . d _{n y} ┘ ^T ∈ ^{n _y}

[0061] are vectors composed of the objective function weights. The vector e is a vector of all ones:

e=[ 11 . . . ] ^T ∈ ^{2n _y} (18)

[0062] and A _u ∈ ^{4n _u xn _u} and b _u ∈ ^{4n _u} are given by 7 $\begin{array}{cc}{A}_u=\left[\begin{array}{c}I\\ -I\\ I\\ -I\end{array}\right]b_u=\left(\begin{array}{c}\stackrel{\_}{u}-u_s\\ -\underset{\_}{u}+u_s\\ N\Delta \stackrel{\_}{u}\\ -N\Delta \underset{\_}{u}\end{array}\right)& \left(19\right)\end{array}$

[0063] while A _y ∈ ^{2n _y xn _y} and b _y ∈ ^{2n _y} are given by: 8 $\begin{array}{cc}{A}_y=\left[\begin{array}{c}I\\ -I\end{array}\right]b_y=\left(\begin{array}{c}\stackrel{\_}{y}-u_s\\ -\underset{\_}{y}+y_s\end{array}\right)& \left(20\right)\end{array}$

[0064] Additionally, since Δy depends linearly upon Δu, the entire problem can be expressed in terms of Δu and ε, reducing the number of decision variables. The resulting LP is cast in standard form and passed to an optimizer. It will be appreciated by those of ordinary skill in the art that the form of A _y , A _u and b _y , b _u given above are specific examples of matrices that can be used in the problem. However, other forms may also be used. For the purposes of the present disclosure, it will be understood that the form of A _y , A _u and b _y , b _u may be more general as dictated by the control purpose at hand in a given application. The present disclosure is intended to address a general form of A _y , A _u and b _y , b _u . It is also to be understood that the slack formulation presented is one of many possible ones currently employed by practitioners in the art. It is believed that the present disclosure is applicable to many other possible forms of implementation.

[0065] The robust LP explicitly accounts for uncertainty in the steady-state gain matrix G of the nominal problem of Equation (15). The mathematical parameterization of possible plants is known as the uncertainty description. It can take many different forms and can be expressed parametrically or statistically. In model identification, if the process noise and model parameters are assumed to be normally distributed or Gaussian variables, the natural uncertainty description is an ellipsoidal bound on the parameters. As illustrated in FIG. 4 , the parameter vector θ=[θ ₁ . . . θ _{n θ} ] ^T ∈ ^{n _θ} is assumed to lie in the set:

θ∈ΘΔ{θ: (θ−θ _c ) ^T V ⁻¹ (θ−θ _c )≦α} (21)

[0066] describing the joint confidence region. The center of the ellipse θ _c (designated with reference numeral 80 in FIG. 4 ) is the mean of the normal distribution and the symmetric matrix V 0 is the covariance (or estimate of the covariance) of the distribution. (The notation “A 0” is meant to be read as: “A is positive definite.”) The matrix V ⁻¹ gives the size and orientation of the ellipsoid. In particular, the square roots of the reciprocals of the eigenvalues of V ⁻¹ are the lengths of the semi-axes of the ellipsoid, and the eigenvectors of V ⁻¹ are the directions of the semi-axes.

[0067] A different form of Equation (21) may in some cases be useful. The confidence region defined by the ellipsoid 80 in Equation (21) and designated with reference numeral 80 in FIG. 4 can be rewritten as:

θ∈ΘΔ{θ _c +α ^½ V ^½ s: ∥s∥≦ 1} (22)

[0068] The term α in the above equations comes from requiring that the parameters be known with a confidence level of γ∈[0, 1].

α=Φ ⁻¹ (γ) (23)

[0069] Here Φ(x) is the univariate normal distribution function with zero mean and unit variance: 9 $\begin{array}{cc}\Phi \left(x\right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^x{}^{-\frac{1}{2}s^2}s& \left(24\right)\end{array}$

[0070] A value of α is always guaranteed to exist for any γ, even for the case of when V is only an estimate of the true covariance (see, e.g., Y. Bard, “Nonlinear Parameter Estimation,” Academic Press, Inc., New York, 1974). For Equation (22) to remain convex, γ≧0.5, or equivalently, α≧0. The uncertainty description given by Equation (22) is more general than Equation (15) as it allows V to be rank deficient.

[0071] When the covariance matrix V does not have full rank, the ellipsoid becomes degenerate, collapsing in specific directions corresponding to parameters which are known exactly. The directions in which the ellipse collapses correspond to the eigenvectors (or semi-axes) with associated eigenvalues that are zero (i.e. those directions for which a nonzero z produces V ^½ z which is zero.) A similar discussion of ellipsoid uncertainty is also given in A. Ben-Tal and A. Nemirovski, “Robust Solutions of Uncertain Linear Programs via Convex Programming,” tech. rep., Technion Institution of Technology, Haifa, Israel, 1994 which gives a more formal treatment of the ellipsoid bounds in the general context of convex programming, and in S. Boyd, C. Crusius, and A. Hansson, “Control Applications of Nonlinear Convex Programming,” Journal of Process Control , (1998), which proposes the use of ellipsoid uncertainty descriptions for use in optimal control. In the steady-state target calculation, the uncertainty parameters θ are the elements of the steady-state gain matrix G. Let g ₁ be the column vector describing the i ^th row of G:

G=└g ₁ g ₂ . . . g _{n y} ┘ ^T (25)

[0072] Assume that outputs are perturbed by independent random variables. No output is correlated to another. In this case a block diagonal covariance matrix can be constructed for the process gains made up of the individual covariance matrices for each row of the gain matrix. In particular, if V ₁ is the covariance matrix corresponding to the i ^th row, then the covariance for the entire matrix is given by 10 $\begin{array}{cc}V=\left[\begin{array}{cccc}{V}_1& & & \\ & V_2& & \\ & & ⋰& \\ & & & V_{n_y}\end{array}\right]& \left(26\right)\end{array}$

[0073] If the outputs are actually cross-correlated, then the off diagonal elements of V will be nonzero. In any case, it can be assumed without loss of generality that g is nominally {tilde over (g)} and has covariance V, yielding:

g∈ξΔ{{tilde over (g)}+α ^½ V ^½ s: ∥s∥≦ 1} (27)

[0074] where ξ defines the ellipsoid for the entire matrix. This is the most natural way to pose the problem as it results in constraint-wise uncertainty. In the nominal LP, the gain matrix G is premultiplied by the (possibly dense) matrix A _y in the output constraint. This results in a linear combination of outputs. The i ^th component of the output constraint is given by 11 $\begin{array}{cc}\sum _ja_{\mathrm{ij}}g_j^T\Delta u\le b_{y_i}& \left(28\right)\end{array}$

[0075] In the general case, ellipsoids of the type above are required to capture the uncertainty in the constraint. If, however, most of the elements of A _y are zero, as is the case for simple bounds on the outputs, it is possible to simplify the elliptic uncertainty description.

[0076] For simple bounds, there is a constant scalar factor β premultiplying the i ^th constraint.

β g _i ^T Δu≦b _{y 1} (29)

[0077] Each component of the constraint depends only upon a single row of the gain matrix. Thus an ellipsoid can be defined for only a single row. Let the gain for the i _th row nominally be {tilde over (g)} _i . The covariance is V _i . The ellipsoid is given by:

g _i ∈ξ _i Δ{g ₁ +α ^½ V _i ^½ s: ∥s∥≦ 1} (30)

[0078] Additionally, there are times when it is desirable to approximate the ellipsoidal uncertainty as either a polytope or simple box.

[0079] If the ellipsoidal constraint can be interpreted in terms of a joint confidence region, then the box constraint can be interpreted in terms of joint confidence intervals (or as approximations to the ellipse). A simple box constraint, designated with reference numeral 82 in FIG. 5 , gives the minimum and maximum bounds of the ellipse. More complicated polytopes may also be used in an attempt to approximate the ellipse. A general polytope, including box constraints, can be expressed as bounds on some linear combination of the elements g _1j of G. The uncertainty description β is then given by:

g∈βΔ{g: A _g g≦b _g } (31)

[0080] Box constraints, when used to approximate elliptic constraints, are more conservative than ellipsoidal constraints and in turn lead to a more conservative control action. The ellipsoid and box constraints described by Equations (27), (30), and (31) each describe a possible uncertainty description U defining realistic values of the gains for use in the robust LP.

[0081] As an illustrative example of the robust LP concept, consider the linear steady-state system:

Δy=GΔu (32)

[0082] with Δu∈ ² and Δy∈ . Let the nominal value of g be: 12 $\begin{array}{cc}\tilde{g}=\left[\begin{array}{c}1.00\\ 0.75\end{array}\right]& \left(33\right)\end{array}$

[0083] Because there is only one output, g=g ₁ and ε=ε ₁ .

[0084] The variance for the problem is: 13 $\begin{array}{cc}V=\left[\begin{array}{cc}0.03& -0.05\\ -0.05& 0.40\end{array}\right]& \left(34\right)\end{array}$

[0085] For the present example, the bias and slack variables are ignored. It is necessary that the matrix elements be known with a probability level γ of 99%, which corresponds to α=2.33. The ellipsoid uncertainty can be approximated as box constraint with the following data: 14 $\begin{array}{cc}{A}_g=\left[\begin{array}{cc}1& 0\\ 0& 1\\ -1& 0\\ 0& -1\end{array}\right]b_g=\left(\begin{array}{c}1.26\\ 1.71\\ -0.74\\ 0.21\end{array}\right)& \left(35\right)\end{array}$

[0086] The ellipsoid uncertainty description and its box approximation are illustrated in FIG. 6 , where the ellipsoid uncertainty description is designated with reference numeral 84 and the box approximation is designated with reference numeral 86 .

[0087] A problem arises when considering output constraints in the nominal problem of Equation (15). The constraints can correspond to critical variables—the temperature limit on a reactor vessel or the maximum allowable feed rate into a unit, for example. Uncertainty or “fuzziness” in the constraint means that even though the nominal value may be binding at the upper or lower limit, the actual value may be outside the constraint region. To see this, the following constraints may be added to the example: 15 $\begin{array}{cc}-4\le \Delta u_1\le 4-3\le \Delta u_2\le 1-2\le \Delta y_1\le 4& \left(36\right)\end{array}$

[0088] The feasible region defined by the constraints corresponds to the interior of the polygon designated with reference numeral 88 in FIG. 7 . This region 88 assumes nominal values for the gains.

[0089] The lower diagonal line, designated with reference numeral 90 in FIG. 7 describes the lower bound on the output Δy ₁ ≧−2. The line 90 is given by:

g ₁₁ Δu ₁ +g ₁₂ Δu ₂ +2=0 (37)

[0090] Uncertainty in g causes there to be uncertainty in the constraint. FIG. 8 shows the nominal constraint and its confidence bounds. Consider the case in which the constraint is binding. Uncertainty could force the true process to lie anywhere within the confidence limits, meaning it could easily lie outside of the feasible set. The equation describing the bounds is given by:

{tilde over (g)} ^T Δu±α ^½ ∥V ^½ Δu∥+ 2=0 (38)

[0091] where ∥·∥ is the standard Euclidean norm. Every output constraint will contain some degree of uncertainty. If the uncertainty from the upper constraint Δy _i ≦4 is included, the plot of the boundary of the new feasible set is as shown in FIG. 9 . The boundary is given by the solid line designated with reference numeral 92 . FIG. 9 shows the original feasible region as the dotted line designated with reference numeral 94 and the box uncertainty approximation of the ellipsoid as the dashed line designated with reference numeral 96 . Box uncertainty will always lead to a more conservative control action than ellipsoid uncertainty. The box uncertainty description contains extraneous gains or plants that are not physically realizable. As the controller is required to consider more and more possible plants, corresponding to an increase in the size of the uncertainty description, and accordingly robustness, the size of the feasible region will decrease, corresponding to a more conservative controller. This is the familiar robustness/performance tradeoff.

[0092] The feasible region can be defined as the set of inputs which produce feasible outputs for any possible value of the gains: 16 $\begin{array}{cc}\Delta u\in \left\{\Delta u:\begin{array}{c}{A}_y\Delta y=A_yG\Delta u\le b_y,\forall G\in U\\ A_u\Delta u\le b_u\end{array}\right\}& \left(39\right)\end{array}$

[0093] where U is one of the uncertainty descriptions discussed above. In addition to gain uncertainty affecting the constraints, it also has an effect on the objective function:

J _s c ^T Δu+d ^T Δy (40)

= c ^T Δu+d ^T GΔu.

[0094] Let the gradient of the objective with respect to the inputs be ƒ:

∇ J _s =ƒ=c+G ^T d (41)

[0095] For the present example, assume c and d are given by:

c=[7.00 5.00] ^T and d=[−6.00].

[0096] Uncertainty causes the value of ƒ to change. ƒ is nominally {tilde over (ƒ)}:

{tilde over (ƒ)}=[1.00 0.50] ^T

[0097] and can take on values anywhere between +E,uns ƒ and {overscore (ƒ)}:

+E,uns ƒ=[−0.56 2.56] ^T {overscore (ƒ)}=[−5.26 6.26] ^T .

[0098] FIG. 10 shows the effect of this uncertainty in the objective function. As the objective function changes, the solution of the LP changes. One solution is to include the gains in the objective, minimizing them in some sense. From a control viewpoint, however, this may not be desirable. An alternative control policy would be to minimize the “best guess” or nominal value of the gain. The goal is to drive the process toward the economic optimum using the best guess for the gain of the process but restrain the inputs to only those that ensure the process will remain feasible for any reasonable gains (those captured by the uncertainty description). If the nominal value of the gain is used, the objective function for the robust LP becomes:

J _s =c ^T Δu+d ^T {tilde over (G)}Δu+e ^T ε (42)

[0099] Although the example is simple, it illustrates the nature of the problem. Model uncertainty forces the problem structure to change. The linear constraints so commonly used in control theory become nonlinear.

[0100] In the foregoing example, it was possible to graphically determine the optimal solution. In practice, there is a need to solve the following optimization problem: 17 $\begin{array}{cc}\underset{\Delta u,\epsilon }{\min }42)c^T\Delta u+d^T\tilde{G}\Delta u+e^T\varepsilon \mathrm{subject}\mathrm{to}& \left(43\right)\\ A_yG\Delta u\le b_y+\varepsilon -A_yb,\forall G\in UA_u\Delta u\le b_u\varepsilon \ge 0& \left(44\right)\end{array}$

[0101] Equations (43) and (44) constitute the robust LP in accordance with the presently disclosed embodiment of the invention. The robust LP is preferably solved at every time-step during the execution of the MPC algorithm. Solution times are preferably on the order of tens of seconds (or less) for it to be used in a real-time setting. Its solution is passed as targets to the dynamic MPC calculation block (block 64 in the exemplary embodiments of FIGS. 2 and 3 ). The problem has been cast in terms of only Δu and ε, having removed the problem dependency on Δy. This has done herein for the sake of simplicity and to illustrate the structure of the problem. Those of ordinary skill in the art will appreciate that the problem can just as well be written so that the outputs appear explicitly as an equality constraint.

[0102] Next, it is desirable to formulate the optimization problem of Equations (43) and (44) as an optimization problem over a second-order cone. To illustrate this, only the simpler problem of the following Equation (45) need be considered: 18 $\begin{array}{cc}\underset{x}{\min }c^Tx\mathrm{subject}\mathrm{to}\mathrm{AGx}\le b,\forall G\in U& \left(45\right)\end{array}$

[0103] which is actually a semi-infinite program. A semi-infinite program generally has the following form: 19 $\begin{array}{cc}\underset{x}{\min }f\left(x\right)h\left(x,\theta \right)\le 0\mathrm{for}\mathrm{all}\theta \in U& \left(46\right)\end{array}$

[0104] with ƒ: ⁿ → and h: ⁿ × ^m → ^p . θ is generally a set of parameters and, for most engineering problems U is closed and convex. The vector x is finite while the set U is infinite. A finite number of variables appear in infinitely many constraints. The problem of Equation (47) occurs in many different fields ranging from robotics to structural design. R. Hettich and K. Kortanek, “Semi-infinite programming: Theory, Methods and Applications,” SIAM Review, 35 (1993), pp. 380-429 provides a review of various areas to which semi-infinite programming has been applied as well as current theory and solution methods.

[0105] The most straightforward way to handle constraints of the type in Equation (47) is to note that:

h ( x,θ )≦0∀θ∈ U (47)

[0106] is equivalent to 20 $\begin{array}{cc}\underset{\theta }{\max }\left(h\left(x,\theta \right),\theta \in U\right)\le 0& \left(48\right)\end{array}$

[0107] If the following definition is used 21 $\begin{array}{cc}H\left(x\right)\underset{\underset{\_}{\_}}{\Delta }\frac{\max }{\theta }h\left(x,\theta \right)& \left(49\right)\end{array}$

θ∈U

[0108] then the semi-infinite program in Equations (46) and (47) becomes: 22 $\begin{array}{cc}\underset{x}{\min }f\left(x\right)H\left(x\right)\le 0& \left(50\right)\end{array}$