Predictive-repetitive control with constraints: From design to implementation

Journal of Process Control 23 (2013) 956–967

Contents lists available at ScienceDirect

Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont

Predictive-repetitive control with constraints: From design to implementation Liuping Wang a,∗ , Chris T. Freeman b , Shan Chai a , Eric Rogers b a b

School of Electrical and Computer Engineering, RMIT University, Melbourne, Victoria 3000, Australia Electronics and Computer Science, University of Southampton, Southampton SO17 1BJ, UK

a r t i c l e

i n f o

Article history: Received 24 January 2012 Received in revised form 19 December 2012 Accepted 31 March 2013 Available online 21 June 2013 Keywords: Repetitive predictive control Frequency sampling filter Model predictive control Robotic arm experiments

a b s t r a c t This paper addresses the design and implementation of predictive-repetitive control systems, including the case when constraints are imposed. The approach is based on a frequency response decomposition technique that detects the dominant frequency components of the reference signals and embeds them in the predictive-repetitive controller. Analysis and design make use of the frequency response of the closedloop system and the choice of the structure of the repetitive controller is considered as a balance between the reduction of tracking errors and minimization of the effects on performance of unwanted elements, such as model uncertainty. The implementation of operational constraints on the amplitudes of the control signals, their increments and on the outputs is considered using an on-line solution constructed using quadratic programming and the identification of active constraints. Experimental results from application to a 2 degree-of-freedom robotic system are used to illustrate the design and implementation procedures developed. © 2013 Elsevier Ltd. All rights reserved.

1. Introduction In many process control applications the controller set-points do not periodically change and control methods such as three term (PID) or model predictive control (MPC) see, for example [1,2] are adequate if any disturbances present are non-periodic, or if the frequency of the disturbance is lower than the bandwidth of the controlled plant. Signals arising in a wide variety of control systems application domains, such as process control, aerospace and robotics, are periodic or can be approximated as such over a large time interval. The general control problem is to force the plant output to track a periodic reference signal and/or reject periodic disturbances. Repetitive control [3] arose from such applications and is based on the idea of using information from previous periods to update the control signal applied on the current one. Previous research [4–6] has developed predictive-repetitive control (PRC) where the frequency components of both the reference and disturbance signals, if present, are analyzed and those selected as dominant embedded in the control law design using the internal model principle [7]. The frequencies to be included are determined by Fourier decomposition of the reference and/or disturbance signals. In this paper, the performance of the closed-loop

∗ Corresponding author. Tel.: +61 3 992 52100; fax: +61 3 992 52007. E-mail addresses: [email protected], [email protected] (L. Wang). 0959-1524/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jprocont.2013.03.012

repetitive control system is examined when such a decomposition is used. Attention is restricted to the reference signal or vector, since the extension to the presence of disturbance signals is a straightforward generalization. The inclusion of operational constraints is also examined, resulting in real-time implementation using a quadratic programming procedure proposed in Ref. [1] that determines the active constraints in real-time by finding the corresponding Lagrange multipliers. Supporting experimental results confirm that it is possible to solve the real-time control problem with a relatively fast sampling rate without generating a high computational load. The remainder of this paper is organized as follows. In Section 2, the PRC algorithm without constraints is summarized; Section 3 discusses how to incorporate constraints in the design, including real-time computation in the presence of both input and output constraints. Section 4 gives supporting experimental results from controlling a robotic arm configured as a two-input and twooutput system. The experimental results also consider the case of conflicting constraints and, in particular, the identification and modification of those active to give a stable controlled system.

2. Predictive-repetitive control without constraints The PRC design considered in this paper is based on a frequency domain decomposition of the supplied reference signal or vector in the single-input, single-output (SISO) and multiple-input, multiple-output (MIMO) cases respectively. Once these are selected

L. Wang et al. / Journal of Process Control 23 (2013) 956–967

they are embedded in the process state-space model in accordance with the internal model principle. These steps are considered in turn below. Consider an SISO example and suppose that a periodic reference signal with period T is uniformly sampled with interval t to give the corresponding discrete sequence r(k), k = 0, 1, . . ., M − 1, where M is defined as T/t and is assumed to be an odd integer. It is also assumed that the periodic reference signal contains no frequency components higher than 1/2t. By Fourier analysis, this discrete periodic signal can be uniquely represented by the inverse discrete Fourier transform (IDFT) as 1 R(ejωi )ejωiik M M−1

r(k) =

(1)

i=0

where R(ejωi ), i = 0, 1, 2, . . ., M − 1, are the frequency components and, for notational simplicity, the fundamental frequency is expressed as ω = 2/M. The z-transform of the periodic signal r(k) is given by Rm (z) =

M−1 

r(k)z −k

(2)

or



(M−1)/2

R(ejlω )F l (z)

(3)

l=−(M−1)/2

where Fl (z) is the lth frequency sampling filter given by F l (z) =

1 1 − z −M 1 = (1 + ejlω z −1 + · · · + ej(M−1)lω z −(M−1) ) M 1 − ejlω z −1 M

Rm (z) ≈ Ra (z) =

n  

Re(R(e

jli ω



))FRli (z) + Im(R(ejli ω ))FIli (z)

+ · · · + d z −

(7)

where D(z) is termed the ‘signal generator polynomial’ for the specified periodic signal. The extension to MIMO examples is immediate and the details are omitted. The controller is to be designed to follow the reference trajectory and hence, by the internal model principle [7], the required D(z) must be included in the denominator of its z transfer-function description. In this paper the requirement is met by adding a vector term to the state dynamics in the plant state-space model as described next. Suppose that the plant to be controlled has mu inputs and my outputs and is represented by the state-space model xm (k + 1) = Am xm (k) + Bm u(k) + (k)

(8)

where xm (k) is an n1 × 1 state vector, u(k) is an mu × 1 input vector, y(k) is an my × 1 output vector and each entry in the n1 × 1 vector (k) is the inverse z-transform of 1/D(z). In practice the action of a zero-order hold means that the input signal u(k) cannot directly affect the output y(k) at the same sample instant and hence the feedthrough term Dm is assumed to be zero for the remainder of this paper. Let q−1 denote the backward shift operator and D(q−1 ) the shift operator interpretation of D(z). Applying D(q−1 ) to xm (k) and u(k) of (8) gives

Also D(q−1 )(k) = 0 and from (8) xs (k + 1) = Am xs (k) + Bm us (k)

2(1 − cos(li ω)z −1 )(1 − z −M ) 1 − 2 cos(li ω)z −1 + z −2 2 sin(li ω)z −1 (1 − z −M ) 1 − 2 cos(li ω)z −1 + z −2

(9)

D(q−1 )y(k + 1) = Cm Am xs (k) + Cm Bm us (k) where the output equation can also be written as y(k + 1) = −d1 y(k) − d2 y(k − 1) − · · · − d−1 y(k −  + 2)

(10)

Introducing the state vector (4)

where Re and Im, respectively, denote the real and imaginary parts of their complex number argument, and l1 , l2 , . . ., ln correspond to the indices of the significant frequency components. The real and imaginary components of the filter are given by

1 M

= 1 + d1 z

−1

− d y(k −  + 1) + Cm Am xs (k) + Cm Bm us (k)

i=1

FIli (z) =

(1 − 2 cos(li ω)z −1 + z −2 )

i=1

1 1 − z −M R(ej0 ) M 1 − z −1 +

1 M

n 

D(z) = (1 − z −1 )

xs (k) = D(q−1 )xm (k), us (k) = D(q−1 )u(k)

If the sampling filters are assumed to have zero initial conditions then the output of each of them in response to a unit impulse input applied at the origin is a periodic signal. If the DFT coefficients corresponding to one or more frequency components is judged to be insignificant then the corresponding filter can be removed from the sum in (3), hence reducing the number used. In the case of n significant frequency components in a periodic signal, the z-transform of Rm (z) is approximated as

FRli (z) =

common denominator of the terms in Ra (z) is given by

y(k) = Cm xm (k) + Dm u(k)

k=0

Rm (z) =

957



x(k) = xsT (k)

yT (k −  + 1)

...

x(k + 1) = Ax(k) + Bus (k)

T

(11)

y(k) = Cx(k) where



A=

The denominator of the z-transform of the zero frequency filter is 1 − z−1 , which is the integrator factor included in discrete-time feedback control systems to, for example, reject constant disturbances or eliminate unknown system bias. Furthermore, each pair of filters for an arbitrary frequency, say li ω, results in the denominator having the factor 1 − 2 cos(li ω)z−1 + z−2 . Hence the least

yT (k − 1)

gives the following augmented state-space model for PRC design

(5)

(6)

yT (k)

Am

0

Cˆ

Ad

⎡ ⎢ ⎢ ⎢ Ad = ⎢ ⎢ ⎢ ⎣





Cm Am

, Cˆ =



0

−d1 I

−d2 I

...

−d−1 I

−d I

I

0

0

...

0

0

I

0

...

0

.. .

..

..

.

..

.. .

0

0

...

I

.

.

0

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

958

L. Wang et al. / Journal of Process Control 23 (2013) 956–967

and 0 and I denote the zero and identity matrices, respectively, of compatible dimensions (my × my ). In addition B

=

C

=

 

(Cm Bm )T

T Bm

0

I

0

...

0 0

... 0



0

0

T

L(m + 1) = ˝L(m) where L(m) = l1 (m)

2.1. Unconstrained PRC design In MPC, modeling the future control trajectory is a critical task. The traditional approach is to embed an integrator in the design and the incremental control trajectory is then directly computed within an optimization window. For PRC, the signal to be optimized is the filtered control signal us (k), and the design can be undertaken by modeling this signal using pulse functions [6]. The main drawback is the requirement to optimize a large number of parameters if fast sampling is required and/or the system has a relatively complex dynamic response [8,9]. Fast sampling is typically required for mechanical and electromechanical systems because the time constants arising in the various sub-components can vary in duration and a smaller sampling interval t is required to capture the effects of the smaller of these. One approach to reduce the number of parameters requiring optimization on-line is to parameterize the future trajectory of the filtered control signal using a set of Laguerre functions, where a scaling factor is used to reflect the time scale of the predictive control system. The use of Laguerre functions in MPC and PRC, including identification based models of the system dynamics, is detailed in, for example [1,8,9] and the following is a summary relevant to the new results in this paper, focusing on the SISO case with the natural MIMO extension noted. In the PRC design, optimization within a moving horizon window is used to find the optimal control signal. Assume, therefore, that the sampling time index k denotes the beginning of the moving horizon window and let the integer m be such that 0 ≤ m ≤ Np , where Np is the length of the moving horizon window, also called the prediction horizon [2]. Then using state-space model (11) m−1 

Am−i−1 Bus (i)

(12)

i=0

where x(k) is the state vector at sampling instant k and x(k + m | k) is predicted state vector at sampling instant k + m, m > 0, given x(k) . The basis of the design is the use of a set of discrete orthonormal functions to describe the filtered future control signal us (m) within a moving horizon window, 0 ≤ m ≤ Np . Assume that N is the number of terms in the expansion and let li (m), 1 ≤ i ≤ N, be a set of Laguerre functions, which are orthonormal. Then

us (m) ≈

N 

(15)



The poles of (11) are the union of those for plant model and that representing the embedded reference signal frequencies.

x(k + m | k) = Am x(k) +

where 0 ≤ a < 1 is the scaling factor. Also the network structure of the z transfer-function can be exploited to show that the set of discrete Laguerre functions satisfies the difference equation

ch lh (m)

(13)

⎡

a

0

..

.

.. .

0⎥ ⎥

..

.

0

0⎥

.. .

..

.

aN−3 ˇ

...

ˇ

a

−a3

···

(−1)N−1 aN−1

ˇ

 h (z) =

1 − a2

1 − a z −1



z −1 − a 1 − a z −1

⎥

⎥ ⎥ ⎥ ⎥ .. . 0⎦

(16)

ˇ = (1 − a2 ) and L(0) =

 

ˇ 1 −a a2

T

Setting a = 0 and ıi (m) = ı(i), where ı(i) is the Dirac delta function that recovers the standard formulation of model predictive control. The extension to the MIMO case follows immediately on considering each input channel and the corresponding column of the matrix B. Continuing with the MIMO case, given the state vector x(k) the prediction of the future state at time m, written x(t + m | k) can be written as x(k + m | k) = Am x(k) + (m)

(17)

where, if Bi is the ith column of B,



= T1

T2

...

Tmu

T

and (m) =

m−1 



T (j) Am−j−1 B1 L1T (j) . . . Bmu Lm u



j=0

Moreover, the number of terms (N) and the scaling factor (a) used in this last construction can be chosen independently for each input. The state vector prediction (12) has a convolution sum that requires the future trajectory of the control vector us (i), 0 ≤ i ≤ Np . The basic idea in Laguerre function based PRC and MPC is to represent us (i) by a set of Laguerre functions together with their unknown coefficients. This is illustrated in the SISO case by us (i) = LT (i)

(18)



where the Laguerre function vector L(i) = l1 (i)



l2 (i)

···

T

lN (i)

T

and the Laguerre coefficient vector = c1 c2 · · · cN . Moreover, N is the dimension of the Laguerre function vector and is also the number of terms used in the approximation. The Laguerre functions are pre-determined in the design once the scaling factor 0 ≤ a < 1 and the number of terms N are chosen [1,8]. Substituting (18) into (12) gives x(k + m | k) = Am x(k) +

m−1 

Am−i−1 BLT (i)

(19)

i=0

and the PRC design problem is find the Laguerre coefficient vector that minimizes the cost function

h−1 (14)

0

and

⎤

...

a

−aN−2 ˇ

lN (m)

...

⎢ ⎢ ˇ ⎢ ⎢ ˝=⎢ ⎢ −aˇ ⎢ .. ⎢ ⎣ .

h=1

where, in general, the coefficients ch are functions of k but the notation here is used for simplicity. In this application, the z transfer-function of the hth Laguerre function, see also Ref. [10], is given by

l2 (m) . . .

T

J=

Np  m=1

xT (k + m | k)Qx(k + m | k) + T RL

(20)

L. Wang et al. / Journal of Process Control 23 (2013) 956–967

where Q is a symmetric positive semi-definite matrix and RL is a block diagonal matrix, where in each block, the diagonal entries are equal to the same positive real number. Further analysis of the cost function, including the extension to the MIMO case [1,8] shows that it can be written as J = T ˝ + 2 T x(k) + E

(21)

959

in, for example [11] and a prescribed degree of stability, that is, the requirement that the controlled system poles lie inside a circle of radius  with || < 1. In the first case, each term in the cost function is multiplied by ˛ to the corresponding power where the choice of ˛ > 1 changes the relative weighting on current and future errors. The prescribed degree of stability case is also straightforward and hence neither of these modifications is detailed in this paper.

where E is a constant, and the matrices ˝ and are given as ˝=

Np 

T

(m)Q(m) + RL ; =

m=1

Np 

2.2. Frequency response analysis (m)QAm .

m=1

Without constraints, the cost function is minimized at sampling instance k when = −˝−1 x(k)

(22)

and by the receding horizon control principle the filtered control signal is us (k) = LT (0)

(23)

where L(0) is the initial Laguerre function vector. As an illustrative example for the two-input two-output case with an identical scaling factor and number of terms



us (k) =

LT (0)

0

0

LT (0)





With no constraints imposed in the MIMO case, substitute (22) into (23) and replace x(k) by the estimated xˆ (k) of (27) and write us (k) in the form us (k) = −Kmpc xˆ (k)

(29)

Taking the z-transform of (27) gives ˆ X(z) = (zI − A + Kob C)−1 BU s (z) + (zI − A + Kob C)−1 Kob Y (z)

(30)

where Us (z) and Y(z) denote the z-transforms of us (k) and y(k), respectively. Hence using (29) the z-transform of the control vector U(z) is U(z) = −

−1 1 {I + Kmpc (zI − A + Kob C)−1 B} D(z)

× Kmpc (zI − A + Kob C)−1 Kob Y (z)

(24)

(31)

or U(z) = F(z)Y(z) where

The control vector to be applied to the system is given by

F(z) =

D(q−1 )u(k) = us (k),

−1 1 {I + Kmpc (zI − A + Kob C)−1 B} Kmpc (zI − A + Kob C)−1 Kob D(z) (32)

or u(k) = us (k) − d1 u(k − 1) − d2 u(k − 2) − · · · − d u(k − ).

(25)

Moreover, state vector x(k) contains the filtered state vector xs (k) and the outputs y(k), y(k − 1), . . . and y(k − ) . At sampling instant k, the tracking errors are e(k)

=

y(k) − r(k) (26)

=

S(z) = (I + Gm (z)F(z))−1

r(k −  + 1) − y(k −  + 1)

For implementation of the control law, these error signals are used instead of the original output entries in x(k) . In many applications all entries in the state vector xm (k) in (8) will not measurable. In such cases an observer is required for control law implementation. Letting xˆ (k) denote the estimated state vector, a full order observer takes the form xˆ (k + 1) = Aˆx(k) + Bus (k) + Kob (y(k) − C xˆ (k))

(27)

where Kob is the observer gain matrix selected such that the eigenvalues of the system matrix A − Kob C are strictly inside the unit circle of the complex plane and at the desired locations. It is also possible to design a reduced order observer for the plant dynamics only. Letting xˆ s (k) be the estimated vector, the reduced order observer is

(34)

where Gm (z) = Cm (zI − Am )−1 Bm is the z-transfer-function matrix of the plant model. Also let ||Gm ||∞ denote the H∞ norm of Gm (z), that is, ||Gm (ejω )||∞ =



∗ (ejω )) (Gm (ejω )Gm

where ( · ) denotes the maximum singular value of its matrix argument and the superscript ∗ denotes the complex conjugate transpose operation. Moreover, by construction, D(z) = 0 at the discrete frequencies 0, ±2/M, ±4 /M, ±6/M where M is the number of samples within one period of the reference signal, and hence ||T(ejω )||∞ = 1 and ||S(ejω )||∞ = 0 at these frequencies. The following analysis investigates how the approximation of the reference signal affects the tracking errors in the closed-loop system. Let E(z) = R(z) − Y (z)

∞

r (ys (k) − Cm xˆ s (k)) xˆ s (k + 1) = Am xˆ s (k) + Bm uˆ s (k) + Kob

(28)

r is selected such that the where the reduced observer gain matrix Kob r C has all eigenvalues closed-loop observer system matrix Am − Kob strictly inside the unit circle of the complex plane. After obtaining xˆ s (k), the full estimated state vector xˆ (k) is constructed using



xˆ (k) = xˆ sT (k)

yT (k)

yT (k − 1)

···

yT (k − )

T

By modifying the cost function, it is straightforward to extend the above analysis to include exponential data weighting ˛, detailed

(33)

or the complementary sensitivity function T (z) = I − S(z)

.. . e(k −  + 1)

The bandwidth of the closed-loop PRC system is characterized by the sensitivity function

(35)

where E(z) = e(k)z −k is the z-transform of the closed-loop k=0 feedback error e(k). Also routine manipulations give E(z) = (I + Gm (z)F(z))−1 R(z) = S(z)R(z)

(36)

The approximation used to construct R(z) means that high frequency components are neglected in the design, but in application these will translate into periodic signals with small magnitudes in / 0, except in the special case where the refe(k). Hence lim e(k) = k→∞

erence signal R(z) contains a limited number of frequencies and all of them are embedded in the PRC.

960

L. Wang et al. / Journal of Process Control 23 (2013) 956–967

The sum of squared errors for one period M can be examined using Parseval’s theorem, see, for example [12], where it is shown that for a periodic signal p(k) with period M

3.1. Formulation of the constraints

 1 2 p(k) = |Pk |2 M

umin ≤ u(k) ≤ umax

M−1

M−1

k=0

k=0

(37)

where the sequence with elements |Pk |2 is the power spectrum of the periodic sequence p(k). Applying (37) in the SISO case gives

  1 |e(k)|2 = |Ek |2 = |Sk |2 |Rk |2 M M−1

M−1

k=0

M−1

k=0

(38)

trace



1 e(k)eT (k) M M−1

k=0

= trace

M−1 

and those on the rate of change on the input as umin ≤ u(k) ≤ umax . Suppose also that an integrator is used in the controller, and the polynomial D(q−1 ) of (25) is expressed as D(q−1 ) = (1 − q−1 )(1 + d1 q−1 + d2 q−2 + · · · + d−1 q−+1 ) Then since

k=0

where the sequence |Ek |2 forms the power spectrum of the periodic error sequence e(k), the sequence |Sk | is formed from the absolute values of the sensitivity function evaluated at the discrete frequencies, k = 0, 2/M, . . ., 2(M − 1)/M and |Rk | is the absolute value of the kth Fourier coefficient of the periodic reference signal. At higher frequencies, if the Fourier coefficients of the reference signal are sufficiently small then the magnitude of the sensitivity function |Sk | becomes relatively large, depending on the controlled system dynamics, where it should be of unity magnitude or slightly larger at some frequencies. Therefore, at both lower and higher frequencies, the approximation used in the design does not lead to significant errors. At median frequencies, the Fourier coefficients of the reference signal could still be significant if the polynomial / 0 at these frequencies and hence errors arise in this freD(z) = quency band. However, the contribution of these to the sum of squared error is determined by |Sk |2 |Rk |2 where the sensitivity function weight |Sk |2 = / 0 should be small due to the extension of the closed-loop frequency bandwidth. In the MIMO case, e(k) is a column vector and the trace, the sum of the squares of its diagonal elements, of the matrix M−1 T e(k)e(k) is used as a quantification of the error, that is, k=0



Suppose that the input amplitude constraints are specified as

us (k) = D(q−1 )u(k) = (1 + d1 q−1 + d2 q−2 + · · · + d−1 q−+1 )u(k) the rate of change of the control signal is u(k) = us (k) −

−1 

L(0) ≤ umax +

 

The wider the bandwidth of a feedback control system the more the measurement noise will be amplified. The attenuation of measurement noise is characterized by the complementary sensitivity function T(z) and hence in the bandwidth where ||S||∞ is zero, ||T||∞ is unity, and the measurement noise in this band will be amplified. Moreover, robust stability of the closed-loop system is also characterized by T(z). In qualitative terms, if the modeling error between the dynamics and its transfer-function description is small when ||T||∞ is large, the closed-loop system will be stable. However, if ||T||∞ is large and the modeling error is also large, the closed-loop repetitive control system could become unstable. The conclusion is that if a larger number of frequencies are to be embedded in the design, more effort should be put into deriving the dynamic model in the higher frequency band.

di u(k − i)

(41)

i=1

−LT (0) ≤ −(umin +

 

di u(k − i))

(42)

di u(k − i)

(43)

i=1

(39)

k=0

(40)

For notational simplicity, the constraints are formulated for SISO systems only in this paper, the extension to the MIMO case is immediate. To impose the constraints at sampling instant k using the moving horizon window, us (k) is related to the parameter vector as given in (23) and us (k + m), m = 2, 3, . . ., can also be expressed using in an identical manner. Hence the input constraints at the sample instant k can be written as

 Sk Rk Rk∗ Sk∗

di u(k − i)

i=1

LT (0) ≤ umax +

−1  i=1

−LT (0) ≤ −(umin +

−1 

di u(ki − i))

(44)

i=1

If transfer-function models are used in the PRC design, the entries in the state vector are estimated using an observer and these do not, in general, correspond to the plant physical variables. Hence state vector constraints will not be imposed as they lack physical meaning when using an observer and it is difficult to quantify their actual bounds. Output constraints are imposed using the predicted output from the augmented model (11) and the one step ahead output is calculated using y(k + 1 | k) = CAx(k) + CBus (k)

3. Constrained PRC In many control systems constraints are imposed to meet physical limits and to ensure safe plant operation. Both input and output constraints are commonly encountered in control applications, where for the former the underlying idea is to express the current and future constraints in terms of the parameter vector us and past states of the controller such that the constrained control problem is solved by minimizing the cost function (20) subject to a set of linear inequality constraints.

(45)

The two-step ahead prediction for the output is y(k + 2 | k) = CA2 x(k) + CABus (k) + CBus (k + 1) = CA2 x(k) + (CABLT (0) + CBLT (1))







(46)

LT b

where L(1) is the Laguerre function for k = 1 and LbT = CABLT (0) + CBLT (1). Moreover, the constraints can be extended to further steps

L. Wang et al. / Journal of Process Control 23 (2013) 956–967

ahead but often actual implementation is based on the first two steps and the rest of the constraints are relaxed in order to reduce the real-time computational load. In this paper, the output constraints are imposed on the second step predicted output using (46). Also e(k + 2 | k) = y(k + 2 | k) − r(k + 2) and the constraint imposed is emin ≤ e(k + 2 | k) ≤ emax Hence in PRC design LbT ≤ emax + r(k + 2) − CA2 x(k)

(47)

−LbT ≤ −emin − r(k + 2) + CA2 x(k)

(48)

≥0

M ≤ 

(49)

⎡ umax

⎤

 

+ di u(k − i) ⎢ ⎥ ⎢ ⎥ i=1 ⎢ ⎥ ⎡ LT (0) ⎤  ⎢ ⎥  ⎢ ⎥ min −(u + d u(k − i)) ⎢ ⎥ i ⎢ −LT (0) ⎥ ⎢ ⎥ ⎢ ⎥ i=1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ −1 T ⎢ L(0) ⎥  ⎢ ⎥ ⎢ ⎥ max + ⎢ ⎥ M=⎢ ;  = u d u(k − i) ⎥ i ⎢ ⎥ ⎢ −LT (0) ⎥ ⎢ ⎥ i=1 ⎢ ⎥ ⎢ ⎥ ⎢ T ⎥ −1 ⎢ ⎥  ⎣ Lb ⎦ ⎢ ⎥ −(umin + di u(k − i)) ⎥ ⎢ T ⎢ ⎥ −Lb ⎢ ⎥ i=1 ⎢ max ⎥ 2 ⎣ e + r(k + 2) − CA x(k) ⎦ −emin

constraints. This formulation also means that at sampling instant k (22) and (25) satisfy all of the constraints. If this is the case, the constrained and unconstrained optimal solutions are identical but if (22) produces a control signal that partly violates the constraints then iterative computation is needed to find the active constraints and hence the solution that minimizes the cost function with these in place. One established approach in the literature is to directly search for the active constraints using optimization of Lagrange multipliers, see, for example [15] and this algorithm is summarized next. Using the Kuhn–Tucker conditions [15,16], the optimization problem is to find the parameter vector by solving the following problem: max min[ T ˝ + 2 T x(k) + 2T (M −  )]

and all the constraints can be expressed in the form

where

961

2

− r(k + 2) + CA x(k)

The PRC control design problem is to minimize the cost function J of (21) subject the linear inequality constraints given by (49) in real-time. Moreover, the M matrix is independent of real-time data and can therefore be computed off-line but  depends on the input and state vectors and must therefore be updated in real-time. Next, the on-line solution of this problem using a quadratic programming algorithm [1,13] is described. 3.2. On-line solution for constrained PRC PRC design poses three major challenges when the solution uses quadratic programming, the first of which is the lack of a guarantee for the existence of a constrained optimal solution. For example, if both the input and output constraints are violated they will be in conflict and the quadratic programming solution becomes infeasible. Secondly, since the solution is computed in real-time and involves an iterative computation, the algorithms used must be sufficiently fast to enable computation within one sampling period. This is one of the main reasons why MPC finds many applications in process control, as discussed in relatively early literature by, for example [14], where the sampling period can be of the order of minutes. In contrast, for electro-mechanical systems such as the robotic system considered in this paper, the sampling period is much less than a second. Thirdly, as for all other engineering applications, there must be safety protection mechanisms in place to ensure stability of the controlled system in the event that the quadratic programming algorithm fails to reach a constrained optimal solution. The quadratic programming algorithm used in this paper has features suitable for real-time computation of this constrained optimal problem in terms of these challenges. Iterative computation in the solution of the quadratic programming problem arises from the inequalities (49) that represent the

(50)



where  is the vector of Lagrange multipliers, ˝ and are the data matrices from the PRC cost function (21) and M and  are given by (49). Also the cost function is minimized with unconstrained when = −˝−1 ( x(k) + M T )

(51)

and, on substituting (51) into (50), the Lagrange multiplier vector  is obtained by solving min[T M˝−1 M T  + 2T (M˝−1 x(k) +  ) + x(ki ) T ˝−1 x(k)] T

≥0

(52) where the maximization problem (50) has been converted to this minimization problem by multiplying the objective function by −1. For notational simplicity, (52) is written as min[T H + 2T K]

(53)

≥0

from this point onwards where H = M˝−1 MT and K = M˝−1 x(k) +  and the final term is independent of  and hence can be neglected. Although this is still a quadratic programming problem that requires an iterative solution, the constraints are much simpler ( ≥ 0) and iterative solutions much easier to find. There are at least two approaches in the literature for solving the optimization problem (53), the first of which is Hildreth’s algorithm [17]. On iteration m + 1 this algorithm solves for the ith entry, denoted i , in the Lagrange multiplier vector  using m+1 = max(0, wim+1 ) i with wim+1 = −

⎡ 1 hii

⎣ki +

i−1 

(54)

hij m+1 + j

j=1

n 

⎤ ⎦ hij m j

(55)

j=i+1

where n is the dimension of the Lagrange multiplier vector, hij is the element in row i and column j of H and ki is the ith element in K. The second approach to solving (53) is to use Parallel Quadratic Programming (PQP) [18]. In this algorithm, a diagonal matrix with non-negative diagonal entries and denoted by diag(rw ), is used to define the matrices H+ and H− as H + = max(H, 0) + diag(rw )

(56)

H − = max(−H, 0) + diag(rw )

(57)

where max (a, b) is applied element-wise. Also define the vectors K+ and K− vectors as K + = max(K, 0)

(58)

962

L. Wang et al. / Journal of Process Control 23 (2013) 956–967

K − = max(−K, 0)

(59)

Suppose that 0 > 0 is an initial guess for the Lagrange multiplier vector. Then on iteration m + 1, the ith entry in the Lagrange multiplier vector is computed as

 m+1 = m i i

Ki− + (H − m )i

 (60)

Ki+ + (H + m )i

where Ki− , (H− m )i and (H+ m )i denote the ith elements of the corresponding vectors. Once convergence is achieved, the resulting  is used to compute the parameter vector and hence the PRC controller using (51) where the constraints corresponding to the non-zero Lagrange multipliers are active and the rest inactive. 3.3. Dealing with conflicting constraints The most commonly encountered cause of conflicting constraints arises when those active are linearly dependent. In common with Hildreth’s algorithm, the PQP algorithm can be shown to converge if the active constraints are not in conflict. However, all quadratic programming algorithms will diverge if the active constraints conflict, see, for example [15] but since the solutions given by the Hildreth and PQP algorithms do not involve matrix inversion this conflict does not immediately interrupt online computation. Moreover, conflicting constraints will result in divergence of the iterative computation of the Lagrange multipliers. As Hildreth’s algorithm is based on the gradient method that uses a basis vector [0, 0, . . ., 1, . . ., 0, 0]T at each iteration m, it attempts to minimize the constrained cost function (53). The advantage of using this algorithm is that when conflicts occur, even though the iterative computation is terminated without convergence, it gives a sub-optimal solution of the Lagrange multiplier vector  for the iteration when this situation occurs. The ability of a quadratic program algorithm to deal with conflict constraints is paramount in the implementation of predictive control. When the active constraints become linearly dependent, that is, conflict, one remedy is to relax one or more of the less critical constraints and this strategy can be implemented by ranking the constraints in terms of their importance. For example, the amplitudes of the control signals used, u(k), are ranked as the most important if they cannot be violated either because of operational limits or because of safety protection of the equipment. Less important are the output constraints via the feedback error, e(k). The changes in the amplitudes of the control signals, u(k), may be less important than the amplitudes themselves, depending on the application. If linear dependence, reflected by the divergence of the Lagrange multipliers, is detected, it is likely that one or more constraints are violated. At this stage, the algorithm should automatically check whether the critical constraints are satisfied and if this is the case the constrained optimal solution will be accepted although the less critical ones are compromised. Conversely, if the critical constraints are violated, the algorithm will eliminate the less critical ones and thereby resolve the conflict. This elimination can be implemented by first setting the Lagrange multiplier vector entries with indices corresponding to the less critical constraints to zero and then form the matrices Mact and  act by choosing the rows in the original M and  of (49) corresponding to the positive Lagrange multipliers. New positive Lagrange multipliers are then found using the closed-form formula

act = −(Mact ˝

−1

T −1 Mact ) (act

+˝

−1

x(k))

(61)

Fig. 1. Picture of the robot arm showing pick and place locations.

and the constrained optimal solution, satisfying the critical constraints in the presence of conflicts, is T = −˝−1 ( x(k) + Mact act )

(62)

4. Experimentally verified design The remainder of this paper deals with the design and experimental implementation of the PRC algorithm of the previous section applied to a laboratory based robot undertaking a pick and place operation that replicates many industrial applications of repetitive and iterative learning control algorithms. Next the required details of the robot are given. Fig. 1 shows an anthropomorphic robot arm undertaking a ‘pick and place’ task in a horizontal plane using two joints. The robots end-effector travels from the ‘pick’ to the ‘place’ location in a straight line using joint reference trajectories that minimize the end-effector acceleration. During the movement, the arm stops at two intermediate points, chosen such that there is a change in the direction of travel along the path after reaching each of them. Having reached the ‘place’ location, the robot repeats the movement in reverse, arriving back at the ‘pick’ location. Positional and velocity control loops have been implemented around each joint to provide baseline performance and the control scheme operates at 20 Hz (t = 0.05 s). Two sets of operational tasks are considered in this paper, both of which are mapped to the coordinates of the x and y axes, respectively, and termed the r1 and r2 reference signals. Fig. 2 shows the r1 and r2 signals and, since t = 0.05 (s), the number of samples within one period is M = 400. Their reconstructions using the first, second and third frequencies, respectively, based on the frequency content data for these signals given in Table 1. Using these Fourier coefficients, the frequency sampling filters models are constructed, see also [6], and their time-domain responses are compared in Fig. 2. Table 2 lists the mean squared Table 1 DFT coefficients of the reference signals for the two reference trajectories. Freq. components

0

1st

2nd

3rd

1 Rm 2 Rm

167.1810 679.8660

−195.6558 −115.2289

61.6788 −146.4782

29.1642 −12.6684

Table 2 Mean squared error of reconstructed signal for the first set of trajectories. Freq. components

0

0–1st

0–2nd

0–3rd

MSE of channel 1 MSE of channel 2

0.5344 0.4394

0.0583 0.2743

0.0110 0.0074

0.0004 0.0054

L. Wang et al. / Journal of Process Control 23 (2013) 956–967 2

963

2.5

1.5 2

1 0.5

1.5

0 −0.5

1

Channel 1 reference

−1

Reconstruction with 0−1st frequency

−1.5

Reconstruction with 0−2nd frequency Reconstruction with 0−3rd frequency

−2 0

2

4

6

8

10

12

14

16

Channel 2 reference Reconstruction with 0−1st frequency Reconstruction with 0−2nd frequency Reconstruction with 0−3rd frequency

0.5 0 18

20

(a) r1 reference signal reconstruction example

0

2

4

6

8

10

12

14

16

18

20

(b) r2 reference signal reconstruction example

Fig. 2. Reference signal reconstruction using a limited number of frequencies.

errors (MSE) for the selection of model orders in the frequency selective filter structure. The control objective is for each output to follow the corresponding reference signal as closely as possible in the presence of operational constraints. To study the effect of the number of frequencies on the bandwidth of the PRC system, the following cases are considered (A) D(z) = (1 − z −1 )

(63)



(B) D(z) = (1 − z −1 ) 1 − 2 cos

 2k  400



(C) D(z) = (1 − z −1 ) 2k=1 1 − 2 cos



(D) D(z) = (1 − z −1 ) 3k=1 1 − 2 cos

z −1 + z −2

 2k  400

 2k  400

 (64)

z −1 + z −2

z −1 + z −2



4.1. Experimental results: PRC without constraints In the experimental results reported in this paper, the parameters for the Laguerre functions are N1 = N2 = 6, and a1 = a2 = 0.95. The weighting matrices are R = I and Q = CT C. The exponential weight is set at ˛ = 1.02, and the prescribed degree of stability is  = 0.95 (see the discussion at the end of Section 2.1). For the observer design, Qob = I and Rob = 40I. Additional disturbance and measurement noise has been added to the robot arm, as shown in Fig. 4.

1

0 (65)

−1 0



50

100

150

200

100

150

200

(66)

1

||S(ejω )||∞

Fig. 3 shows the Bode plot of for the four cases over 0 ≤ ω ≤ . At the lower frequencies, the Fourier coefficients are significant, see Table 1, the magnitude of the sensitivity function |Sk | is sufficiently small, see Fig. 3, and is zero at the frequencies where D(z) = 0.

0

−1 0

50

Time (sec)

0.1 0

0

10

−0.1 0

−5

10

50

100 Time (sec)

150

200

50

100

150

200

0.1 −10

10

0 −15

10

−0.1 0 −4

10

−2

10

0

10

Time (sec)

Frequency (rad) Fig. 3. Plot of the H∞ norm of the sensitivity function S(z) for Case A (*); Case B (+); Case C (o); Cased D (square).

Fig. 4. From top to bottom: input disturbances as voltage fluctuation on u1 ; input disturbances as voltage fluctuation on u2 ; output measurement noise on y1 ;output measurement noise on y2 .

964

L. Wang et al. / Journal of Process Control 23 (2013) 956–967

reference signal

reference signal 1.5 Output y (t) (rad) 1

Output y1(t) (rad)

1.5 1 0.5 0

0

−0.5

−0.5 0

20

40

60

80

100 Time (s)

120

140

160

180

0

200

20

40

60

80

100 Time (s)

120

140

160

180

200

0

20

40

60

80

100 Time (s)

120

140

160

180

200

0

20

40

60

80

100 Time (s)

120

140

160

180

200

0

20

40

60

80

100 Time (s)

120

140

160

180

200

2.5

2.5

2

Output y (t) (rad)

2 1.5

1.5

2

Output y (t) (rad) 2

1 0.5

1 0.5

1 0.5 0

0 0

20

40

60

80

100 Time (s)

120

140

160

180

200

Error e1(t) (rad)

0.6 0.4 Error e1(t) (rad)

0.2 0

−0.2

0.5

0

−0.4 −0.5

0

20

40

60

80

120

140

160

180

200 1

Error e (t) (rad)

0

0.5

2

2

Error e (t) (rad)

0.5

100 Time (s)

0

−0.5

−0.5 0

20

40

60

80

100 Time (s)

120

140

160

180

200

Fig. 5. Experimental results for Case A. The output positions are given in the top two plots, and the output position errors in the bottom two.

Fig. 6. Experimental results for Case B. The output positions are given in the top two plots, and the output position errors in the bottom two.

Consider Case A, which is equivalent to the use of a predictive controller with integrators and embedding no periodic components of the reference signals into the design. The experimental results are shown in Fig. 5 and it is evident that there are steady-state errors between the corresponding reference and the output position signals. In Case B, Fig. 3, confirms that the bandwidth of the PRC system is widened with the sensitivity function equal to zero at 2/400 (rad). Fig. 6 shows the experimental results for this case, which confirm that both y1 and y2 follow their reference signals r1 and r2 without steady-state errors, as confirmed by the error plots shown in the same figure. As the bandwidth of the feedback control system is further increased by including another periodic mode in the design, that is, Case C defined by (65) the closed-loop performance degrades, as confirmed by Fig. 7. The oscillations that appear in the closed-loop responses are likely to be caused by the unmodeled high frequency dynamics in the robot arm. In addition, the measurement noise is amplified. Finally, when all three periodic modes are included in the design model, that is, Case D defined by (66) the system becomes unstable. At this stage, a single set of performance tuning

parameters have not been found to stabilize the closed-loop system, and the instability of this system is caused by the neglected high frequency dynamics of the robot arm. 4.2. Experimental results: PRC with constraints In the second set of experimental results given in this paper, operational constraints are added to the unconstrained design of the previous sub-section. The first case considered is when constraints are placed on the input amplitudes only and the second when constraints are placed the input amplitudes, their increments and the output errors. Case I. The constraints are −0.3 ≤ u1 (k) ≤ 1.6, 0 ≤ u2 (k) ≤ 2.5

(67)

The experimental response plots of the resulting control system are shown in Fig. 8 and confirm that constraints are satisfied with all of them active. The evidence for this observation is that the control signals at various stages reached their maximum and minimum values. This confirms that the closed-loop performance does not deteriorate due to the active input operational constraints.

L. Wang et al. / Journal of Process Control 23 (2013) 956–967

2

965

2

y1

Output y (t) (rad)

1.5 1

r1

1

0.5

1

y1min

0

0

−0.5 −1 −1.5 0

20

40

60

80

100 Time (s)

120

140

160

180

200

−1 3

0

10

20

30

40

50

60 y2

3

r2

2

2

y2min

2

Output y (t) (rad)

y1max

reference signal

1

1 0

−1 0

20

40

60

80

100 Time (s)

120

140

160

180

200

0

10

20

30 Time (sec)

40

50

60

2

1

1

0

u1

Error e1(t) (rad)

y2max

0

−1

0 0

20

40

60

80

100 Time (s)

120

140

160

180

200

−1

0 3

10

20

30

40

50

60

10

20

30 Time (sec)

40

50

60

2

1 0

u2

Error e2(t) (rad)

2

−1 −2

1 0

0

20

40

60

80

100 Time (s)

120

140

160

180

200

0

Fig. 7. Experimental results for Case C. The output positions are given in the top two plots, and the output position errors in the bottom two.

0.2

Case II. The input amplitudes, the difference in the input signals and the output signals are constrained as follows:

Δu1

0.1 0 −0.1

−0.3 ≤ u1 (k) ≤ 1.6

(68)

0 ≤ u2 (k) ≤ 2.5

(69)

−0.5 ≤ (r1 (k) − y1 (k)) ≤ 0.5

(70)

−0.5 ≤ (r2 (k) − y2 (k)) ≤ 0.5

(71)

−0.12 ≤ u1 (k) ≤ 0.12

(72)

−0.2

−0.12 ≤ u2 (k) ≤ 0.12

(73)

−0.4

−0.2

It is much more difficult to stabilize the robot arm with all the constraints imposed compared with the previous case where only those on the input amplitude were imposed. The computational demand is not a significant issue with respect to the sampling interval of 0.05 s. With both input and output constraints imposed, the system is unstable when both become active and start conflicting. This situation was caused by the external disturbances added to

0

10

20

30

40

50

60

0

10

20

30 Time (sec)

40

50

60

0.4

Δu2

0.2 0

Fig. 8. Case I results: top two plots outputs and errors, middle two plots inputs with constraints marked, bottom two plots differences in inputs with the constraints marked as horizontal dashed lines.

966

L. Wang et al. / Journal of Process Control 23 (2013) 956–967

2

y1

1

r1 y1min

0

y1max

−1 0

10

20

30

40

50

60

3

y2

2

r2 y2min

1

y2max

0 0

10

20

30 Time (sec)

40

50

60

2

u1

1 0 −1

0

10

20

30

40

50

60

output constraints represented by (70) and (71) are less important and, for this robot arm, are ranked second. The constraints on the incremental control signals u1 (k) and u2 (k) (72) and (73) are used to reduce voltage spikes and are allowed to be momentarily violated over a very short interval (violation over longer periods in many applications would trip the protection equipment, halting the experiment). For this application, they are ranked as third in importance and in application the least important constraints in a conflict situation will be relaxed first with the others remaining active. Using the solution of the quadratic programming problem using the active Lagrange multipliers, the experimental implementation, after detecting conflicting constraints, sets the Lagrange multipliers that correspond to the constraints to be relaxed equal to zero, according to the rank of importance. This resolves the problem of conflicting constraints and leads to the convergence of the active Lagrange multipliers. Fig. 9 shows the corresponding experimental results and it is evident that the controlled system is stable. Both input amplitude and output constraints are satisfied but the differences in the control signals, u, have violated their constraints on the occasions when conflicts have occurred and they were chosen to be relaxed.

3

5. Conclusions

u2

2 1 0

0

10

20

30 Time (sec)

40

50

60

40

50

60

0.2

Δu1

0.1 0 −0.1 −0.2 0.2

0

10

20

30

This paper has undertaken an experimental evaluation of constrained PRC applied to a robot arm, where the real-time optimization is performed by solving the active Lagrange multipliers. In the case of conflicting constraints, the algorithm relaxes one or more of them to ensure the convergence of the quadratic programming computation. A representative set of experimental results have been given and discussed for the cases of input amplitude constraints alone and for constraints on input amplitudes, their differences from sample-to-sample and the output signals. These demonstrate that the PRC design considered can be experimentally implemented and, as required, could be considered for industrial deployment.

Δu2

0.1

References

0 −0.1 −0.2

0

10

20

30 Time (sec)

40

50

60

Fig. 9. Case II results: top two plots outputs and errors, middle two plots inputs with constraints marked, bottom two plots differences in inputs with constraints marked as horizontal dash lines.

the system, see Fig. 4, which could not be measured and predicted in advance and therefore could not be avoided through off-line computation. Instead, the real-time optimization algorithm should include a mechanism to detect the presence of conflicting constraints and relax one or more of them to ensure that the closed-loop feedback control system is stable in the extreme case. Therefore, for this application the constraints need to be classified into hard and soft before the experimental verification commences. It is also useful to establish a relative importance rank for the constraints in the event of conflict. In this experiment the manipulated variables are demand voltages sent to the robotic arm. These have magnitude limits which cannot be violated due to risk of damaging the equipment and hence must be treated as hard constraints. From this point of view, the constraints represented by (68) and (69) are ranked first in importance. The

[1] L. Wang, Model Predictive Control System Design and Implementation Using MATLAB, 1st ed., Springer, London, 2009. [2] J.M. Maciejowski, Predictive Control with Constraints, Pearson Education Limited, 2002. [3] S. Hara, Y. Yamamoto, T. Omata, M. Nakano, Repetitive control system: a new type servo system for periodic exogenous signals, IEEE Transactions on Automatic Control 33 (7) (1988) 659–668. [4] L. Wang, P. Gawthrop, D.H. Owens, E. Rogers, Switched linear predictive controllers for periodic exogenous signals, International Journal of Control 83 (4) (2010) 848–861. [5] L. Wang, S. Chai, E. Rogers, Predictive repetitive control based on frequency decomposition, in: American Control Conference, Baltimore, USA, 2010, pp. 4277–4282. [6] L. Wang, S. Chai, E. Rogers, C.T. Freeman, Multivariable repetitive-predictive controllers using frequency decomposition, IEEE Transactions on Control Systems Technology 20 (6) (2012) 1597–1604. [7] B.A. Francis, W.M. Wonham, The internal model principle for linear multivariable regulators, Applied Mathematics and Optimization 2 (2) (1975) 170–194. [8] L. Wang, Discrete model predictive control design using Laguerre functions, Journal of Process Control 14 (2004) 131–142. [9] G. Valencia-Palomo, J.A. Rossiter, C.N. Jones, R. Gondhalekar, B. Khan, Alternative parameterisation for predictive control: how and why? in: American Control Conference, San Francisco, 2011, pp. 5175–5180. [10] B. Wahlberg, System identification using Laguerre models IEEE Transactions on Automatic Control 36 (1991) 551–562. [11] L. Wang, Use of exponential data weighting in model predictive control, in: 40th IEEE Conference on Decision and Control, Orlando, FL, 2001, pp. 5175– 5180. [12] D.G. Manolakis, V.K. Ingle, S.M. Kogon, Statistical and Adaptive Signal Processing, McGraw-Hill. [13] Y. Wang, S. Boyd, Fast model predictive control using online optimization, IEEE Transactions on Control Systems Technology 18 (2) (2010) 267–278.

L. Wang et al. / Journal of Process Control 23 (2013) 956–967 [14] S.J. Qin, T.A. Badgwell, An overview of industrial model predictive control technology., in: Preprints of Chemical Process Control-CPC V, California, 1996. [15] D.G. Luenberger, Linear and Nonlinear Programming, 2nd ed., Addison-Wesley Publishing Company, 1984. [16] S. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004.

967

[17] C. Hildret, A quadratic programming procedure, Naval Research Logistics Quarterly 4 (1957) 79–85. [18] M. Brand, V. Shilpiekandula, C. Yao, S.A. Bortoff, A parallel quadratic programming algorithm for model predictive control, in: 18th IFAC World Congress, vol. 18, part 1, Milan, Italy, 2011, pp. 1031–1039.