Consider a dynamical system described by the state and output equations

where x(t)∈Rnx  is the system state, u(t)∈Rnu is the control input and y(t)∈Rny is the controlled output at time t, and f_c(., .) and h(.) are possibly nonlinear algebraic functions. The MPC approach aims to drive the output y(t) of the system to a chosen reference value r by computing, and applying the control input u(t) such that it minimizes a suitable cost function over a finite time horizon. A typical cost function considered for minimization in real time is a quadratic cost function of the form.

where Q_c is a positive semidefinite matrix, R_c is a positive definite matrix, l_c(.) is a positive semidefinite function defining the terminal cost, t₀ is the initial time, t_f is the final time, and x¯,u¯ is an operating point of the system such that

The state x(t) and the input u(t) are usually required to satisfy certain constraints which may be expressed as x(t)∈X and ut∈U where X and U are convex sets.

Since the MPC approach envisages a real-time solution of the optimization problem, in practice, a DT approximation of the problem is considered choosing a suitable sampling period, T and the problem is solved every time step. Typically, the system dynamics in (1) are approximated by a DT description

where k∈Z denotes the discrete time step 0, 1, 2, …, x(k) = x(kT) and y(k) = y(kT) denote the values of the state and the output at the beginning of the time step k, u(k) gives the control input to be applied during time step k, and f(., .) is the equivalent or approximate state-update function. Further, assuming a uniform sampling, the cost function in (2) is approximated by a predicted cost function of the form

where Q and R are appropriately defined positive (semi) definite matrices, l(.) is a positive semidefinite function used to specify the terminal cost and N is the horizon length. The subscript k + i|k attached to a variable indicates the value of the variable at time k + i predicted at time k.

The finite-horizon MPC optimization problem can be stated as

subject to the constraints

The basic steps in an MPC scheme based on problem (5) are summarized in Figure 1.

Figure 1: 
Steps in a basic finite-horizon nonlinear MPC scheme

If the system is a linear time-invariant system, the DT state equation in (3) can be obtained exactly, typically assuming a zero-order hold at the input, i.e., assuming u(t) = u(k) for kT ≤ t < (k+1)T. Then, if the terminal cost l(.) is a convex quadratic function, problem (5) can be formulated as a convex quadratic programming (QP) problem and solved efficiently. For a general nonlinear system, however, problem (5) is often a challenging problem, particularly for systems with fast dynamics. On the one hand, an exact DT model of system dynamics is usually not easy to obtain, and approximate models lead to inaccurate multi-step predictions of states. On the other hand, even with a relatively more accurate DT model, the state update equation is typically nonlinear and results in a non-convex optimization problem which is difficult to solve in real time except in cases with small horizon lengths.

A magnetic levitation system comprises an object of magnetic material which is suspended in the air under the effect of an electromagnetic field created by a current-carrying coil. A schematic of the system is shown in Figure 2 where the vertical displacement of the object, y, increases as it moves downwards from the top. By varying the current in the coil, the vertical position of the object can be varied. Given the desired vertical position, r, the controller and the current drive generate the necessary current in the coil, based on the measured actual position (and speed) of the object obtained from the sensor (which is typically an optical sensor), to move the object to the desired position.

Figure 2: 
Schematic diagram of a magnetic levitation system

Nominally, the motion of the object is described by the differential equation

where m is the mass of the object, κ is the relevant coefficient of viscous friction and g is the acceleration due to gravity. F_m (y, i) is the position-dependent force on the magnetic object due to the current in the coil and is given by

where L and a are relevant constants (see, e.g., [33, pp. 31-32]).

Defining x₁ = y and x2=y˙ as state variables and u = i as the input, the corresponding state-space description of the system can be expressed as

The state variables and the control input are typically required to satisfy the constraints

A steady-state operating point (x¯,u¯) of (6) has the form

To formulate a reference-tracking MPC problem of the form (5) for system (6) under the constraints in (7), a DT approximation of the state equation in (6) is needed. One can either use the DT equivalent of the state equation in (6) linearized about the operat-ing point associated with the reference point or a nonlinear DT approximation of the state equation in (6) for a small sampling period by using the simple Euler or some Runge-Kutta (RK) method.

While the use of a linearized model results in a QP-based MPC optimization problem, the approximation involved in the linearization leads to prediction errors for initial states not close to the operating point. The use of some approximate discretization, on the other hand, results in a nonlinear DT system model and hence a non-convex MPC optimization problem with a large horizon length, which can be computationally demanding. In the sequel, we explore an alternative approach to formulate the pre-dictive control problem for system (6) with a linear DT approxi-mation of the continuous-time state equation based on a trans-formation of the input variable.

Define a new variable

and rewrite system (6) as

where Ac=010-κ/m,Bc=0-La/2m,cc=0g and C=[1   0]. The state equation in (8) is linear in the transformed input variable v. This transformation, however, has implications on the discretization of the model, on the form of the constraints and the cost function, and on the implementation of the control inputs.

Assuming a constant v throughout the sampling interval T, the DT counterpart of (8) is given by

where A=eAcT,B=∫0TeAcθdθBc and c=∫0TeAcθdθcc.

The assumption of a constant v throughout the sampling period in this DT model means that the actual control input u to be applied to the system is a function of the value of y(t) within the period. This has implications on the control implementation which we shall discuss in the sequel.

Lemma 1: Given a constant v, let the control input u(t) in (6) for any time t∈[t0,t0+T) be defined as

where q=12κ2mg-Lav,q1=mκq-x2(t0) and q2=a+x1t0-q1. Then, the state of system (6) at time t₀+T is

Proof: This result can be easily seen with some algebraic exercise. With v defined as u² / (a + x₁)², the second component of the state equation in (6) can be written as

Assuming a constant v(t) = v for t from t₀ to t₀ + T, we have,

The first component of the state derivative in (6) can be ex-pressed as

The value of x₁(t) for any t∈[t0,t0+T] is given by

Hence, with u(t) defined as in Equation (10), v(t) is constant and equal to v and hence the value of x(t₀ + T) directly follows from the DT equivalent (9) of model (8).

Lemma 1 shows that if the control input u(t) given by Equation (10) is implemented for t∈[t0,t0+T) with v = v(k) for each t₀ = kT, k = 0, 1, 2, …, then the controlled system is exactly described by the DT LTI system model (9). So, under this assumption, we can use model (9) to predict future states of the system.

Using (9) as the model of the system, we can easily express the sequence of the future states of the system at the sampling instants. Defining the vectors

we can express x(k) as (see, e.g. [2])

where

and c = B(c). Here, I_nx denotes an identity matrix of size n_x.

Constraints on state variables in (7a) can be imposed at the sampling instants and are collectively expressed as a linear con-straint

where M_x and M_v are appropriately defined matrices. The constraint on the input variable u, expressed in terms of v, is, however, more complicated. Like y(t), for a constant v(k), u(t) varies within the sampling period and we can impose the constraint on it at the sampling instants, i.e.,

for i = 0, 1, …, N - 1. While the left-hand side constraints in (12) are equivalent to

the right-hand side constraints are nonlinear and non-convex. Define constant matrices

where blkdiag(Δ₁, Δ₂) denotes a block diagonal matrix with diagonal blocks Δ₁ and Δ₂, and ⊗ denotes the Kronecker product. Further, defining a new matrix variable

the constraints in (12) can be equivalently expressed as a pair of linear matrix inequality (LMI) constraints on v(k) and V(k) as

where diag(v) for a vector v represents the diagonal matrix with the elements of v along its diagonal and diag(Δ) for a square matrix Δ represents a diagonal matrix containing the diagonal elements of Δ. Further, BB=B(B) and a = a1_N, where 1_N denotes an N-dimensional vector with all 1’s. Matrix elements marked * in (15) complete the symmetric matrices.

Next, we consider a cost function of the form (4) with Q > 0, a small R ≥ 0 and l(ϕ) = ϕ^T P ϕ where P ≥ 0. This cost function can be written as

where x¯=1N+1⊗y¯0, v¯=1N⊗v¯ with v¯=2mg/La, Q=IN⊗Q00P and R=IN⊗R. Further, omitting the constant terms, we can write J(k) as

The MPC optimization problem can now be stated as the minimization of the cost J(k) in Equation (16) under constraints (11), (14) and (15). Condition (14), however, is a nonconvex equality constraint and is computationally difficult to handle. For the particular form of the cost function that we have, a potential alternative is to convexify it by replacing it with its semidefinite relaxation (see, e.g., [34]) given below:

For the equality constraint (14) to be satisfied, the relaxed condition (17) should hold tightly.

Lemma 2: If the minimization of the cost function J(k) in (16) under constraints (11), (15) and (17) results in an optimal solution such that the relaxed constraint (17) is tight (i.e., condition (14) holds), then then nonlinear constrains in (12) are satisfied by the predicted input and state sequences at all sampling instants.

Proof: It is easy to see that the first set conditions in (12) can be equivalently written as

This can be expressed in vector form as

where ∘ denotes the Hadamard (elementwise) product and the squaring on the left side is to be understood in the elementwise sense. With the assumption that equality (14) holds true, this condition can be written as

which is equivalent to the first LMI in (15). Similarly, the second set of conditions in (12) can be equivalently expressed as the second LMI in (15). Hence, the lemma follows.

Remark 1: In the problem that we have, matrix A and vector c have all positive elements and BB has negative elements. The upper limit constraint on the control input becomes active when x₁(k) is larger than x¯1 and the difference between them is significant. Under this situation, the minimization of cost J(k) in Equation (16) typically results in a solution for which the relaxed condition (17) is tight.

In order to ensure that the output of the controlled system un-der the control scheme based on the minimization of the MPC cost in Equation (16) accurately tracks the desired position, we introduce a terminal constraint in the MPC problem. The terminal constraint set is based on a carefully chosen terminal controller. An offline-optimized dynamic terminal controller such as the one discussed in [11] can offer a larger terminal set which reduces the online computational burden. In this work, however, we consider the standard static feedback terminal control law

where v¯=u¯2/a+x¯12=2mg/La, and K is the optimal state-feedback gain that minimizes an infinite-horizon linear quadratic cost, typically obtained as

where P is the solution to the well-known algebraic Riccati equation and also satisfies the Lyapunov equation

so that the infinite horizon terminal cost is represented by l(ϕ) = ϕ^T P ϕ.

The corresponding terminal constraint set that is computed should be invariant and constraint-admissible for the terminal controlled system, i.e., for a state of the terminal controlled system starting in this set, the system constraints should be satisfied at all future times. Since the control law also depends on x¯, i.e., on x¯1, for terminal set computations, we consider an augmented state z=(x,x¯1), which follows the dynamics

where we have used the fact that Bv¯=-c. Further, since the upper bound on the transformed input variable v, viz., um2/a+x12, is dependent on the state of the system in a nonlinear way, we consider a piecewise affine approximation of this upper bound for a number of segments of x1∈[0,ym] as illustrated in Figure 3. We then compute a constraint-admissible invariant set for the state of system (20) for x¯1 in each of these segments while allowing x_1 to extend beyond the segment under the affine approximation of the upper bound for v for this segment. Let Si,i=1,..,mz be the nonoverlapping segments considered in [0,ym] and, for each segment i, let Zi be the constraint-admissible, invariant set for z(k), satisfying

Figure 3: 
Piece-wise affine approximation of x₁-dependent upper bound on input variable v

obtained for x¯1∈Si and x1 in some S^i⊃Si under the linearized constraint for this section. Note that with the linearized state-dependent bound, the constraints on v(k) can be expressed as linear constraints for each segment and the corresponding invariant set Zi can be obtained using an iterative algorithm (see, e.g., [35]).

The MPC optimization problem including the terminal con-straint based on the simplified prediction strategy for the magnet-ic levitation system can be stated as

with J(k) defined in Equation (16), subject to the constraints (11), (15), (17) and

where Zi is such that x¯1∈Si. The MPC scheme based on problem (22) is summarized in Figure 4.

Figure 4: 
Steps in an MPC scheme based on the simplified prediction strategy

Proposition 3: With the off-line computations carried out as outlined in the scheme of Figure 4, if the MPC optimization problem (22) is feasible for the initial state x(0) at time k = 0, and providing that the constraint (17) holds tightly at every step k ≥ 0, the online steps outlined in the scheme of Figure 4 ensure that (a) the problem remains feasible at all time steps and the system constraints are satisfied at all discrete time instants, and (b) providing that (A, Q^1/2) is observable, the state-input pair of the controlled system is asymptotically driven to the operating point (x¯,u¯) and hence the output y(k) asymptotically tracks the given reference r.

Proof: If problem (22) is feasible with the initial state x(k) at time k, and an optimal solution is obtained with constraint (17) holding tightly, system constraints and the terminal constraint are satisfied by the predicted states and inputs at times k + i, i = 0, .., N - 1 and by the predicted terminal state x(k + N|k). If the input u(t) is applied to the system at times t ∈ [kT, kT + T) as mentioned in the scheme of Figure 4, the state x(k + 1) at time k + 1 is equal to the predicted state x(k + i|k) and the new optimization problem at time k + 1 is feasible with the input variables vk+ik+1=v*k+ik, i=1,..,N-1 and vk+Nk+1=v¯+Kxk+Nk-x¯, and with the predicted terminal state satisfying (xk+N+1k+1,x¯)∈Zi because of the invariance property (20) of the set Zi. Hence, if problem (22) is feasible at time k = 0, part (a) of the proposition follows.

Next, to see part (b) of the proposition, note that since P satisfies the Lyapunov equation (19b), with the optimal state and input sequences predicted at time k, which remain feasible at time k + 1, we have

and hence the cost monotonicity condition

holds. This condition can be shown to ensure part (b) of the proposition using standard arguments (see, e.g., [1], [36]).

The online steps in the MPC scheme of Figure 4 assume that the state x(k) can be measured exactly at every sampling instant and that the time varying control input u(t) can be implemented according to Equation (10). In practice, however, the state may not be measured fully (e.g., only the position of the magnetic object may be measured) or there may be errors in the measurement and the control input may be implemented in a more simplified form, using a constant or a linear approximation. The errors due to these factors can be incorporated in the system description (9) by adding a process disturbance term w(k) in the state equation and a measurement noise term η(k) in the output equation.

Under certain assumptions on the process and measurement disturbances, the state may be estimated using a Kalman filter. Even in the absence of reliable stochastic models of disturbances, the method can be used together with some learning or tuning mechanism that estimates the disturbance covariances (e.g., [37]). However, we consider a receding horizon state estimation using weighted batch least squares [38] to estimate the state. Consider a backward state propagation model based on (9), viz.,

where A~=A-1. Using this model, the stacked vector of states at the current and the past N_e time steps, i.e.,

can be written as

where A~ and B~ are defined just like A and B but with matrix A replaced by A~,c~=B~A~c and

So, the stacked vector of outputs based on x~k is given by

where C=INe+1⨂C. Defining the vector of measured outputs

the error vector is ek=yk-y~k. A weighted least squares estimate of x(k) that minimizes the weighted sum of errors e(k)TWe(k), with a weight matrix W, is given by

where M=A~TCTW CA~-1A~TCTW.

The presence of disturbances in the system dynamics model complicates the design of any controller, and some additional mechanism is usually adopted (e.g., [39]) to ensure robustness. In the MPC approach, a robust mechanism may be followed to guarantee the satisfaction of constraints in the presence of dis-turbances. However, assuming noises of small magnitudes, and relying on the inherent robustness of the MPC approach, we use the prediction model based on the nominal dynamics (9) as dis-cussed in Section 3.2.

For the implementation of the control input, we consider two alternative approximations to the time varying input u(t) for time t ∈ [kT, kT + T) – a linear approximation of (10) given by

where β0=q2+q1e-κT2m1-κT2m and   β1=q-q1κme-κT2m, and a constant approximation based on the time-averaged value of y(t) within the sampling period, given by

where the values of the constant quantities Λ₀, Λ_x and Λ_v are obtained from the definitions of q, q₁ and q₂. When the constant approximation (25b) is used, the two constraints in (15) can be replaced by a single matrix inequality constraint with the left-hand-side matrix having the off-diagonal block equal to

where Lx=IN⨂Λx, Lv=IN⨂Λv and L0=1NΛ0. This version of the constraint applies to the constant control input that is actually applied to the system during each sampling period.

Remark 2: In the analysis here, we have not incorporated the time delay in the implementation of the control action every time step due to the time spent in solving the MPC optimization problem. The effect of this time delay can be accounted for by either solv-ing the optimization problem for the predicted future state of the system or by incorporating the delay in the discrete-time system matrices as mentioned in [2].

For the simulations in this section, we consider the values of the system parameters mentioned on the left side of Table 1.

The limits considered in the constraints in (7) and the cost matrices used in defining the cost J(k) in Equation (16) are mentioned on the right side of Table 1. The terminal cost matrix P, the feedback gain K and the corresponding terminal invariant sets are determined as mentioned in Section 3.3. We consider only two segments S1=[0, 0.05] and S2=[0.05, 0.1] in [0, 0.1] and find the invariant sets Z1 and Z2 for system (20) with x¯1 in S1 and S2. The set Z1 obtained for x1∈[0, 0.06] and T=0.04s is shown in Figure 5.

Figure 5: 
Terminal invariant set Z1  obtained for system (20) with x1∈[0, 0.06] and x¯1∈[0, 0.05]

We first assess the accuracy of the simplified multi-step predictions based on model (9) when linear and constant approximations (25a, b) of the optimal control inputs are applied to the plant. We consider the MPC optimization (22) for a target position r=x¯1=0.0025 with x0=(0.095,0) and compare the predicted values of the state components with their actual values obtained when the inputs predicted at time k = 0 are implemented using the nonlinear relation (10) and using its linear (25a) and constant (25b) approximations within each sampling period. The actual values obtained using u(t) as per (10) obviously coincide with the predicted values. The comparison results shown in Figure 6 for a sampling period T = 0.04s indicate that the actual values in the rest two cases are close to the predicted values through the full prediction horizon. With a smaller sampling period, the actual values remain further closer to the predicted values in both the cases.

Figure 6: 
Position and velocity components of the state based on an exact nonlinear implementation of the control input in (10), a linear approximation of the input (25a) and a constant approximation of the input (25b) for T = 0.04s.

Figure 7 shows the magnitudes of the largest errors, w_{max, i}, i = 1, 2 in the update of the state components in (9) for various values of the sampling interval T when the control input is applied based on the approximations of (25a) and (25b). Note that, in this example, while the error in the update of x_1 is smaller when the linear approximation is used, the error in the update of x₂ is smaller when the constant approximation is used.

Figure 7: 
Largest modelling error component magnitudes in the state equation with linear and constant approximations of u(t).

We also assess how the accuracy of the prediction approach explored in this paper compares with the accuracy of the predictions using the usual linearized model. We consider r = 0.0025 and linearize the nonlinear model about the operating point with x¯=(0.0025,0. We discretize the linearized model for T = 0.04s and consider a finite horizon MPC problem with N = 20 for an initial state (0.095, 0). In Figure 8, predicted trajectories of the state components based on the optimal input sequence obtained at time k = 0 are compared with the actual trajectories obtained by applying this input sequence to the nonlinear system. It can be observed that the predictions are significantly off from actual values beyond a few steps. Note that while a receding horizon control based on these inaccurate predictions may still give a satisfactory control performance, there is no guarantee of stability and accurate set-point tracking with an MPC scheme based on these predictions.

Figure 8: 
Predicted optimal state component trajectories obtained at k = 0 based on the linearized model about a target operating point, and the actual trajectories obtained with the same input sequence

Next, we evaluate the setpoint tracking performance with the overall MPC scheme outlined in Figure 4. Considering x(0) = (0, 0), we first choose r = 0.095 for up to 0.8s and then change the reference value to r = 0.0025. In the simulations, we consider, T = 0.04s and N = 10. Assuming that only the position is measured, we consider a zero-mean Gaussian sensor noise of standard deviation σ_η = 10^-4 in the simulations, and to estimate the state, use a receding horizon state estimator with N_e = 5 and x^0=0,0. A diagonal weight matrix W is chosen with diagonal elements Wii=1/em,i2 where e_{m, i} is an estimated upper bound of the magnitude of e(k - i|k) based partly on the observation in Figure 7 for the constant approximation (25b) and on an assumption that |η(k)| ≤ 2σ_η most of the times.

Figure 9 shows the actual and estimated trajectories of the state components under MPC when the control input u(t) is implemented using the constant approximation (25b). With a small sensor noise, the estimated state components are sufficiently close to the actual state components. Note that the position x₁ reaches and remains close to the relevant set point before any change in the reference input. Also, the constraints on the state components are always satisfied.

Figure 9: 
Actual and estimated position and velocity trajectories of the controlled system under the MPC scheme of Figure 4 for x(0) = (0, 0) and r = 0.095 for up to time 0.8s and r = 0.0025 from time 0.8s. Initial estimated state x^0=(0,0).

Figure 10 shows the actual control inputs applied to the sys-tem. The piecewise constant input trajectory of Figure 10 (left) pertains to the case when (25b) is used to approximate the time-varying control input whereas the piecewise affine input trajecto-ry of Figure 10 (right) is obtained when (25a) is used.

Note that when (25b) is to be used, the two constraints in (15) are replaced by a single constraint as mentioned in Section 3.4. The input trajectories in both the cases satisfy the input con-straints.

Figure 10: 
Control inputs applied to the system under the MPC scheme of Figure 4 when the constant approximation (25b) of (10) is used (left), and when the linear approximation (25a) of (10) is used (right).

The reformulated MPC optimization problem is convex and is solved as an SDP problem. The computational requirement for an SDP problem is typically higher than that for a QP problem which arises in the case of a linearized system with linear con-straints. Nevertheless, the convex reformulation provides a com-putationally efficient alternative to solve the original nonlinear MPC problem which is non-convex. Simulations in this section are carried out using the CSDP solver in MATLAB 2017b in a machine with 1.8GHz Intel i7 processor and 16 GB RAM. The mean computation times needed by the solver at each step in the simulations discussed in the previous subsection (with the input approximation of (25b)) are listed for three horizon lengths in Table 2.

These computation times are small for a small N but grow with increasing values of N. Since the sampling interval considered for the problem is T = 0.04s, if a large horizon length is needed for the problem, computational delay may need to be handled either with a more efficient solver for the specific structure of the problem or accounted for in the DT model of the system as mentioned earlier in Remark 2.