Convergence of the discrete-time nonlinear model predictive control with successive time-varying linearization along predicted trajectories/Diskretinio laiko nuspejamosios kontroles netiesinio modelio ir laiko kitimo teisiskumo konvergencija nuspejamosiose trajektorijose.

Introduction

Significant attention has been given in the last decade to nonlinear model predictive control (NMPC). Currently there are some successful methods [1-4], and applications in nonlinear systems also with finite time horizon and reviews concerning NMPC methods [5-7].

[FIGURE 1 OMITTED]

The nonlinear system described by the discrete-time nonlinear state space model can be rearranged into the so-called state and control dependent linear form [8, 9]. The non-linear behaviour of the system is included in the state and control dependent matrices. If the trajectory prediction for the system may be obtained within the algorithm then one can pretend that the future behaviour is known during the prediction horizon [3]. Such a system can be treated as a linear time-varying (LTV) one. Most often the algorithm has the following common steps shown in fig. 1 [2, 10]. The control can be computed using arbitrary method for LTV systems. Also the technique presented in [3, 4] uses similar idea to [2], but with a different model representation and an optimisation technique.

The main aim of this paper is to analyse convergence of the NMPC successive model linearization method along predicted state and input trajectories. Particularly stopping and necessary convergence condition are discussed.

The algorithm from Fig. 1 refer only to one time step computation. Usually it is employed with receding horizon, where the algorithm must be repeated for successive time steps [k.sub.0] = [k.sub.0] + 1.

Model description

General discrete-time (DT), time-varying nonlinear model is assumed in the following form

x(k +1) = f(x(k),u(k),k). (1)

The nonlinear system can be transformed into following discrete-time, time-varying state-dependent form

x(k+1)= A(x(k),u(k),k)x(k)+B(x(k),u(k),k)u(k), (2)

where state and input dependent matrices are calculated for given initial condition [x.sub.0] and control trajectory u(k) at each time instant.

Then, using the past trajectory, matrices A(k) = A(x(k),u(k),k),B(k) = B(x(k),u(k),k) may be calculated for the subsequent points of the trajectory and the nonlinear system (1) is approximated by the LTV model with matrices A(k), B(k).

DT-LTV system is given in the state space form

x(k +1) = A(k)x(k) + B(k)u(k), (3)

where k = [k.sub.0], [k.sub.0] + 1,..., [k.sub.0] + N -1, A(k)[member of] [R.sup.nxn], B(k)[member of] [R.sup.nxm] and N is the prediction horizon.

It can be equivalently defined using evolution operators or, in the considered finite horizon case, also by following block matrix operators L,B:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

where [[PHI].sup.k.sub.i] = A(k)A(k-1)...A(i). For vectors x, u we use the following block vector notation, i.e.

[??] = [[[x.sup.T]([k.sub.0] + 1)...[x.sup.T]([k.sub.0] + N)].sup.T]. (6)

It follows that the mathematical model can be rewritten in the final form as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

We assume that at each time instant the system can be analyzed as starting from time sample equal to zero with a current initial condition [x.sub.0] = x([k.sub.0]) up to N steps into the future (prediction horizon).

The operator LB is a compact and Hilbert-Schmidt one from [l.sub.2] into [l.sub.2] and boundedly maps signals u(k)[member of]L = [l.sub.2][[k.sub.0],[k.sub.0] + N-1] into signals x [member of] X.

For simulation purposes we employ cost function in following form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)

where P [member of] [R.sup.(nN)x(nN)], Q [member of] [R.sup.(mN)x(mN)] are diagonal weighting operators, constructed with weighting matrices P (k) [member of] [R.sup.nxn], k = 1,2,...,N, Q(k)[member of][R.sup.mxm], k = 0,1,..., N - 1

Convergence of the algorithm

Definition 1. The algorithm from Fig. 1 is convergent if there exists a limiting control sequence [u.sub.opt] such that for any arbitrarily small positive number [epsilon]>0, there is a large integer I such that for all i>I, [parallel][u.sub.(i) - [u.sub.opt][parallel] [less than or equal to] [epsilon]. The algorithm that is not convergent is said to be divergent.

The algorithm converges both for local or global optimal solutions. Divergent algorithm cannot satisfy a stopping condition usually given by following absolute tolerance condition: [parallel][[u].sub.(i)] - [u.sub.(i][parallel] [less than or equal to] [epsilon] for arbitrarily small [epsilon].

Definition 2. Let the state trajectory deviation norm from the reference trajectory for any given time horizon be given by [parallel]x(i) - x[ref][parallel] for the i-th iteration of the algorithm and [x.sub.(i+1)] - [x.sub.ref][parallel] for the next i+1 iteration. The state error rate is denoted by R and defined as the following relation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Theorem 1. The state trajectory in the consecutive iteration of the algorithm from Fig. 1, calculated for any nonlinear system transformed into the state-dependent LTV form under the assumption [parallel][L.sub.(i)][[DELTA].sub.A][parallel] < 1 can be calculated from the following equation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (10)

denotes the difference system operator and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] denotes respectively: the state, the input trajectory, the input operator and the system operator in the i iteration of the algorithm, [x.sub.ref] denotes the reference trajectory.

Proof. The proof follows from similar results for perturbed systems [10] with the difference that the term [[??].sub.A(i)] represents deviation of the linearized time-varying system matrix. The deviation corresponds to corrections of the state and input trajectories which are applied in the consecutive iteration of the NMPC algorithm. Previously, this methodology was used for uncertain systems for which [[??].subA(i)] represented model uncertainty. The system in the i+1 iteration can be treated as the perturbed system from the i iteration, and the system in the i+1 can be written in following form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)

or equivalently

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

To derive the trajectory [[??].sub.(i+1)], the term (I - [L.sub.(i)][[??].sub.A(i)] has to be invertible. Under the sufficient condition that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] becomes invertible. Calculating the left side inverse of above equation leads to eq. (9).

Corollary 1. The state error rate of the algorithm from Fig. 1 can be evaluated from following expression

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (13)

Corollary 2. The state error rate of the algorithm for [[??].sub.ref] = 0 can be evaluated from expression:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

Proof follows directly from the triangle inequality for norms.

Definition 3. The input signal trajectory [??] is called an optimal control trajectory ([[??].sub.opt]) if and only if the cost function (8) attains its global minimum at [??] = [[??].sub.opt]. The optimal state trajectory [[??].sub.opt] is calculated from eq. (1) for given initial conditions. This definition can be applied also for J close-to-global minimum, not only to the exact global minimum.

Let [??] be vector valued function, which transform given linear time-varying system realization in i-th iteration of the algorithm into optimal control trajectory in the i+1-th iteration of the algorithm, such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

Corollary 3. The algorithm from Fig. 1 with transformation g is divergent if the optimal solution [[??].sub.opt] is not the stationary point for vector field [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] defined for fixed f,[x.sub.0] and g.

Proof. The proof follows from definition 2 and the stopping condition. If the condition is not satisfied the algorithm cannot stop even when the optimal control has been found (e.g. by a chance).

The necessary condition for the convergence of the algorithm from Fig. 1 can be described as follows:

Theorem 2. The algorithm from Fig. 1 is convergent in terms of Definition 1 if and only if following condition is held

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

where g is the transformation from given linear time-varying system realization into optimal control trajectory, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are linear time-varying system operators for fixed initial condition [x.sub.0], the optimal control trajectory [u.sub.opt] and the nonlinear function f.

This condition is satisfied for most physical systems but not for all (see numerical example). The NMPC algorithm with a given transformation from eq. (1) into eq. (2) is convergent in terms of definition 1 and theorem 2 if the optimal control [[??].sub.opt] follows directly from the time-varying system operators [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] linearized at [[??].sub.opt].

Proof. Let us assume that there exist optimal input trajectory [[??].sub.opt] for k = 0...N -1 such that

[x.sub.opt](k +1) = A([x.sub.opt](k),[u.sub.opt](k),k)[x.sub.opt](k)+ +B([x.sub.opt](k),[u.sub.opt](k),k)[u.sub.opt](k), (18)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)

The algorithm is convergent to the optimal control trajectory [u.sub.opt] if and only if the following limit exists [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] or equivalently

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

where [u.sub.(i)] is the input trajectory in the i-th iteration of the algorithm.

Let us define iterative control differences vector field in following way

[V.sub.gfX0]([[??].sub.(i)])= [[??].sub.(i+1)]([u.sub.(i)])-[[??].sub.(i)], (21)

where [[??].sub.(i+1)]([[??].sub.(i)]) is the input trajectory in the i+1 iteration of the algorithm given by eq. (16) where initial condition x0, nonlinear function f and transformation g are fixed. Substituting (21) into (20) results in following limit

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (22)

It means that [[??].sub.opt] must be stationary point of the field [V.sub.gfx0]. Taking account eqs. (16) and (21) it may be written

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (23)

It is equivalent to (17) and finishes the proof.

The algorithm is monotonically convergent if R [less than or equal to] 1. The closer the coefficient to zero, the faster the rate of convergence. However, when approaching the optimal solution we get R [right arrow] 1. For linear systems we have R = 1 since the solution is calculated in the first iteration of the algorithm. For values R > 1 the algorithm can be divergent. In some cases R can oscillate above and below 1. In such case the algorithm is convergent if for the consecutive iterations [absolute value of R(i)-1] is decreasing function for i [greater than or equal to] I, where I is a large finite integer. The convergence to the optimal solution J [right arrow] [J.sub.opt] is always connected with approaching zero by nonlinearity differences [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

What could be done to ensure that R is near 1? There are two conditions which have to be satisfied. The operator product L[[DELTA].sub.A] should approach zero or equivalently [parallel][??][[DELTA].sub.A][parallel][right vector]0. From the definition of operator [L] it follows that [parallel][L][parallel] [greater than or equal to] 1. On the other hand the norm of [[DELTA].sub.A] can be arbitrary small. Its actual value depends on the nonlinearity of the system i.e. the nonlinearity degree and the method which decomposes nonlinearity into matrices A and B (step 2 of the algorithm from Fig. 1). For linear systems the matrix A does not depend on the input and state and hence the norm is equal to zero. The assumption A=0 results in [[[DELTA].sub.[??]] = 0. However, it also increases the second difference [L.sub.(i)]([[B].sub.(i+1)][[u].sub.(i+1)]- [[B].sub.(i)][[u].sub.(i)], especially, for small values of input u . In most cases it results in the divergence of the algorithm, similarly as does the assumption B=0. Thus the difference operator norms [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] should be possibly small, and balanced.

Numerical example--convergence necessary condition

It is assumed that the control is calculated iteratively using cost function (8) with [x.sub.ref] = 0, from the formula

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (24)

The algorithm from Fig. 1 is combined with above equation. The convergence of the algorithm is considered for the following dynamical nonlinear discrete-time system

[x.sub.k+1] = [x.sup.2.sub.k] + [u.sup.3.sub.k], [x.sub.0] = 8. (25)

The system can be transformed into the state space dependent form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (26)

The cost function is assumed in following form J = [x.sup.T]x (Q=0). The control in equation (25) is in the 3rd power and the state in 2nd so the optimal control trajectory can be found by hand. It has the following dead-beat form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (27)

Matrices A, B and the system operators for the time horizon N=2 are as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (28)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (29)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (30)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (31)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (32)

The new control [u.sub.n] have following form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (33)

The new control cannot be computed, because the term [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is not invertible on the optimal control trajectory. The necessary condition from theorem 2 is not satisfied and the algorithm is divergent.

Conclusions

Methods proposed in the paper concerns the transformation method from a general nonlinear form into the state space dependent form. The suitability of the chosen transformation method follows from the necessary condition for convergence, what can be deduced from theorem 2 and also from theorem 1, concerning the uniform convergence.

From a practical point of view, the chosen method is suitable if: assumption of theorem 2 is satisfied - the method is not divergent and nonlinearities are decomposed into two additive terms - state and input dependent matrices of the state space dependent form so as to the norms of difference operators of the system [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are possibly small.

References

[1.] Qin S. J., Badgwell T. A survey of industrial model predictive control technology // Contr. Eng. Pract.--Elsevier, 2003. - No. 7(11).--P. 733-764.

[2.] Kouvaritakis B., Cannon M., Rossiter J. A. Non-linear model based predictive control // Int. J. Control.--Taylor & Francis, 1999.--No. 10(72).--P. 919-928.

[3.] Dutka, A., Ordys A. The Optimal Non-linear Generalised Predictive Control by the Time-Varying Approx // Proc. MMAR 2004.--Miedzyzdroje, 2004.--P. 299-304.

[4.] Ordys A. W., Grimble M. J. Predictive control design for systems with the state dependent non-linearities // Proc. SIAM Conf. on Control and its Appl.--San Diego 2001, California.

[5.] Camacho E. F., Bordons C. Model Predictive Control. Springer, 2004.

[6.] Morari M., Lee J. Model predictive control: Past, present and future // Comput. Chem. Engi.--Elsevier, 1999.--No. 4/5(23).--P. 667-682.

[7.] Mayne D. Q., Rawlings J. B., Rao C. V., Scokaert P. O. M. Constrained model predictive control: stability and optimality // Automatica.--Elsevier, 2000.--No. 6(26).--P. 789-814.

[8.] Huang Y., Lu W. M. Nonlinear Optimal Control: Alternatives to Hamilton-Jacobi Equation // Proceedings of the 35th IEEE Conf. on Dec. and Control.--Kobe, 1996.--P. 3942-3947.

[9.] Orlowski P. Nonlinear into state and input dependent form model decomposition--App. to discrete-time MPC with successive time-varying linearization along predicted trajectories // Proc. of the 7th Int. Conf. on Inf. in Control, Aut. and Robotics.--Funchal, 2010.--Vol. 3.--P. 87-92.

[10.] Orlowski P. Convergence of the optimal non-linear GPC method with iterative state-dependent, linear time-varying approximation // Int. Workshop on Ass. and Fut. Dir. of NMPC.--Freudenstadt-Lauterbad, 2005.--P. 491-497.

[11.] Orlowski P. Fractional indexes impulse responses approximation for discrete-time Weyl Symbol computation // Electronics and Electrical Engineering.--Kaunas: Technologija, 2010.--No. 8(104).--P. 9-12.