- Email: [email protected]
Universal algorithm for the solution of the regulator equation
Copyright © IFAC Nonlinear Control Systems, Stuttgart, Germany, 2004
ELSEVIER
PUBLICATIONS www.elsevier.coml1ocnfldifac
UNIVERSAL ALGORITHM FOR THE SOLUTION OF THE REGULATOR EQUATION Branislav Reh:lk .,1 Sergej Celikovsky··,2
• Czech Technical University, Faculty ofElectrical Engineering, Department afControl Engineering, Technickti 2, 166 27 Praha 6, Czech Republic. E-mail: [email protected] .. 1nstitute ofInfonnation Theory and Automation, Academy of Sciences ofthe Czech Republic. P. O. Box 18. 182 08 Prague, Czech Republic. E-mail: [email protected].
Abstract: The computational aspects of the output regulation problem are considered. The standard solution of the output regulation problem uses the explicit solution of the so-called regulator equation being a highly complex PDE combined with algebraic equations. At this moment, the methods to solve regulator equation are, nevertheless, restricted to special classes of systems only. A universal algorithm to solve the regulator equation is proposed here. It is based on a special regularization procedure combin~d ~ith finite elements method and an optimization of the error induced by the regulanzauon. The approach is demonstrated on an example of the two-cart with an inverted pendulum system. Copyright © 2004 IFAC Keywords: Nonlinear control systems, output regulation, partial differential equations.
autonomous exosystem with possible unknown initial states.
I. INTRODUCTION A central problem in control theory and applications is to design a control law to achieve asymptotic tracking with disturbance rejection in a nonlinear system. When a class of reference inputs and disturbances are generated by an autonomous differential equation, this problem is called nonlinear output regulation problem, or alternatively, nonlinear servomechanism problem (lsidori and Byrnes, 1990). The corresponding autonomous differential equation is usually called as the exogenous system. In the sequel, the above setting will be referred to as the output regulation problem. In other words, the output regulation problem treats a possible unknown reference signal and/or disturbances generated by the known neutrally stable
For linear systems the classical output regulation was extensively studied in (Franwon and Murray, 1976; Francis, 1977). For nonlinear systems, the problem was first studied in (Hepburn and Wonham, 1981), and solutions to the output regulation of nonlinear systems have been presented in (Huang and Rugh, 1990; Isidori and Byrnes, 1990) using "full information" which includes the measurements of exogenous signals as well as of the system state. The necessary and sufficient conditions for the existence of a local full information solution of the classical output regulation problem are given in (Isidori and Byrnes, 1990; Huang and Rugh, 1990); they basically mean that the linearized system is stabilizable and there exists a certain invariant manifold. The classical output regulation via error feedback has been solved in (Byrnes et al., 1997; Isidori, 1995) by application of system immersion technique. The plant uncertainty parametrized by
Supported by the Research program J04/98:212300013. Supported by the Gram Agency of the Czech Republic through the research grant No. 10210210709. 1
2
159
unknown constant parameters is treated as a special case of exogenous signals and the solution, extended from the error feedback regulation, is referred to as the structurally stable regulation in (Byrnes et al., 1997).
2. PRELIMINARY FACTS First, ~')me basic facts on nonlinear output regulation problem are recalled (Isidori and Byrnes, 1990; Huang, 2000; Huang, 1995), the notation from (Huang, 2(00) is adopted throughout this paper. Consider the plant
Beginning with (Isidori and Byrnes, 1990), the basic approach to all kind of output regulation problems is to use the explicit solution to the so-called regulator equation which is to be obtained off-line, based on the model of the plant and exosystem only. Nevertheless, regulator equation is a partial differential equation (PDE) combined with algebraic restrictions. Its solvability is a highly complex and questionable issue. From the PDE theory point of view, this is due to its first order structure, not fitting into usual PDE frameworks and, moreover, singularity in its coefficients.
x
= f(x(t)) + g(x(t))u(t) + p(x(t))v(t))
y(t) =
(1)
h(x(t)),
where sufficient smoothness of the vector fields f, g, p and row functions of h is assumed. Further, x(t) E Rn is the state, u(t) E Rm is its input, y(t) E RP its output and w(t) E RV is the so-called exogenous signal. This signal is generated by the so-called exosystem which is supposed to be known and for linear, i.e. for a known (jj x jj )-matrix S and a known (v x jj) matrix Q the exosystem is given by
Existing results solving regulator equation use geometrical approach and special structure of nonlinear systems, for which the solution is trivial. The most simple situation is when overall system (i.e. plant to be controlled plus exosystem) is hyperbolically minimum phase, see already (lsidori and Byrnes, 1990), while (Huang, 2003) shows that for a fairly general class of systems regulator equation may be reduced to its PDE part only. Algorithm to solve this PDE for both minimum and nonrninimum phase cases are suggested in (Huang, 2000; Huang, 1995) based on an undetermined power series technique. Nevertheless, it requires laborious symbolic computations, that are difficult to be implemented as an universal computeraided algorithm applicable automatically to any system from a reasonable class.
iJ = Sv, w(t) = Qv.
(2)
Moreover, the exosystem is assumed to be neutrally stable. This implies the system is Lyapunov stable in both the forward and the backward time direction. Thereby, exogeneous signal is used to describe both reference to be tracked and undesired disturbance to be rejected. This leads to the output regulation problem, which may be tackled by various kind of feedback compensators. The solution to the so-called full information output regulation problem consists in finding the feedback compensator u = a(x, v) such that
The aim of this paper is to provide a universal algorithm for solving regulator equation without referring to some particular structure of the given system. It is based on the regularization procedure that adds new artificial small second order term (Laplacian with small parameter) combined with finite elements methods. The error induced by such a regularization is treated via optimizing certain error functional. The results presented in this paper are encouraging, the proposed algorithm is compared with (Huang, 2000; Huang, 1995) algorithm based on series expansion. It provid es slightly better quality in computing the error zeroing manifold, moreover, this is "automatical" algorithm applicable to any system directly. The theoretical convergence conditions are discussed as well and may be further refined.
(1) the equilibrium x = 0 of the controlled system
is exponentially stable if no exogenous signal is present (2) there exists a neighborhood U C Rn+IJ. of (0, 0) such that for each initial condition (x(O), v(O)) holds limt-++oo e(t) = 0, where the so-called error e(t) is defined as e(t) = h(x(t)) - w(t). The so-called error feedback output regulation problem uses a dynamical feedback compensator using measurements of the error e(t) = h(x(t)) - w(t) only. In both cases, the solution to the following regulator equation is used
8x(v) --a;;-Sv = f(x(v)) + g(x(v))u(v) h(x(v)) = Qv
The paper is organized as follows. Some preliminary facts are presented in the next section, while section 3 presents our algorithm in detail. Section 4 presents illustrative example of the two-cart with an inverted pendulum system. Some concluding remarks are collected in the final section.
(3)
with the condition x(O) = O. More precisely, suppose (x(v), u(v)) is the solution of the regulator equation (3). Then the manifold (x(v), v) is the output zeroing manifold of the augmented plant (1-2) having output h(x) - Qv when with input u = u(v). Based on this property, the feedback compensator solving the full information output regulation problem may be constructed as follows
u = a(x, v) = L(x - x(v))
160
+ u(v),
(4)
bounded domain n c Rn with Lipschitz boundary. The differential equation
L being matrix of gains stabilizing the linear approximation of I. Interested reader may find detailed exposition in (Isidori, 1995), including the error feedback case relying on the solution of regulator equation as well.
8u -a6u + b. 8v has a solution u
3. ALGORITHM FOR SOLVING TIIE REGULATOR EQUATION
V
82
i=1
8v 1
= I: J'i::2
C -
~ div b ~ K > O.
(6)
Corollary: Assume the equation (5) is scalar with k 1 > 0, k 2 > O. Further assume the exosystem (2) is linear and neutrally stable. Then the equation (5) has solution. Proof. The Lyapunov stability means that all eigenvalues have zero real part and none of them is multiple. Hence Trace S = O. Since (using the notation from the lemma above) a = k 1 , b = Sv, C = k 2 we have div b = Trace S = O. Thus the condition (6) is satisfied. D A similar result may be obtained even if the Laplacian 6 is not added. Nevertheless, adding the smallparameter Laplacian is useful since our computational experience shows that the finite-element solution of this augmented equation is "smoother" - exhibits less oscillations on the same mesh without destroying accuracy. Application of the finite-element method for solving the equation is proposed. An extensive treatment of the convergence analysis of numerical schemes for equations like the system (5) is carried out in (Roos, et aI., 1996) and is omitted here. Also, for the sake of brevity investigations of integrability properties of the solution are not presented here. Nonetheless, these properties are crucial for existence of the error functional (see below) as well as for differentiating it. Instead it is assumed all these manipulations are possible.
It consists of two steps - first the control u(v) is chosen and the partial differential equation is solved. Then the control u is adjusted minimizing certain error functional which reflects algebraic part. Then this procedure is repeated until a sufficiently precise approximation of the solution of the system is found. To be more specific, some basic facts on PDE solvability are recalled. Unfortunately, neither existence of its solution nor its uniqueness can be guaranteed for PDE 3. To achieve their solvability their regularization is proposed:
8x(v) -k 1 6x(v) + a:;;-Sv = -k 2x(v) + j(x(v)) + g(x(v))u(v)
.. j, 6
if
a > 0,
The regulator equation (3) has been usually solved using decompositions of the involved functions into the Taylor series in v (Huang, 2000; Huang, 1995). In this paper, such an approach is being replaced by the direct numerical solution of this partial differential equation. As a matter of fact, the system (3) consists of n PDE equations and p algebraic ones in n + m unknown variables x, u. First, without going into details, one may reduce these equation to the case PDE-algebraic system having the same number of equations and variables. The corresponding procedure is similar as the one of (Huang, 2003), nevertheless, not all algebraic equations are to be eliminated in our paper, so that the algorithm of our paper may treat wider class of systems.
+ cu =
(5)
If the error caused by the additional terms in (5) is too large it might be compensated by using the following procedure. Let the symbol r(v) denote the term
with constant vectors k 1 , k 2 with positive (or at least nonnegative) elements.
8x(v) r(v) = a:;;-Sv - (J(x(v))
+ g(x(v))u(v)).
The next issue for (3) is lack of boundary conditions, the only condition posed on the function x is x(O) = O. As it is not possible to solve numerically (5) on the whole space RV, one has to restrict the solution to a sufficiently large bounded domain n c RV and homogenous Dirichlet boundary conditions is used on its boundary. Such a choice of the boundary conditions does not seem to influence the solution significantly if the distance of the trajectory of the exosystem from the boundary of the domain n is large enough and the constants ki are not too large. Thus the condition x(O) = 0 is also satisfied with a fairly good accuracy. Nevertheless, if this condition is unacceptably violated adjusting the boundary conditions is unavoidable.
by adjusting the control u. For the domain holds ~ n. It need not be equal but the set should contain all possible trajectories of the exosystem. The area around the border of the domain n where the boundary condition influence the solution of the differential equation strongly can be out of interest in the error computation.
Theorem 1 «Roos, et aI., 1996) - Lemma 1.18) Let the differential equation in (5) be scalar on a
A formula for the derivatives of the functional J with respect to the control u can be found. The directional
The goal now is to minimize the error functional J =
~
J
(r(v)f r(v)+
(7)
n
(h(x(v) - Qvf(h(x(v)) - Qv)dv
n
161
n n
derivative of the function x(v) in the direction
d dt
X2 X3 X4 Xs
oD x(v)
Sv - kI'~:~.D
x2 h(X2,X3,X.. ,Xs) X4 fs(X2' x3, X4, xs) x6 K -(Xl - xs) M
Xl
~v
x6
(8)
+
(11)
o 1
where Dj, Dg, Dh denote the Jacobi matrices of the functions j, g, h, respectively.
with
h
(9)
can be implemented. Here,
II
u
o o
oD x(v) D
J
o
- COSX3 l(M + m(sinx3)2)
The fact that these derivatives can be computed will enable to utilize some gradient methods for the optimization. For the derivative of the function r holds:
D
+ m(sin X3)2
M
The zero dynamics is represented by the last four equations. The equations describing the motion of the pendulum build up a hyperbolic zero dynamics while the equations of motion of the added cart are nonhyperbolic zero dynamics.
(10)
We seek a control for the system of equations (11) such that the state Xl tracks the exogenous signal v(t) = Asinwt. This signal is generated by the system
Further steps are analogous as in the previous case. The function \11 is adjusted now instead of the function
u.
iJ
= Sv, S = ( _Ow
~) , Q =
( 1 0) .
The equation of the center manifold (3) attain for our example the form:
4. EXAMPLE: lWO-CART WITH AN INVERTED PENDULUM SYSTEM
V2 0X3 (V) _ V1 0X3 (v) =X4(V) OVl OV2 OX4(V) OX4(V) w2 V2--!:l- - V I - - ! : l - = -l Vl COSX3(V) UVI UV2
The method is demonstrated on the example introduced in (Devasia, 1999) which was extensively studied in (Huang, 2(00). The system consists of two carts which are elastically connected. An inverted pendulum is placed on one cart. The input is the force F acting upon the cart carrying the pendulum. The output is the position of this cart, denoted by Xl. The model of the plant is
TSinX3(V) V2 oxs(v) _ OVl
162
Vl
oxs(v) = X6(V) OV2
(12)
+ (13) (14) (15)
The symbol ~i[cpj] stands for the derivative of the function Xi, i = 5,6 with respect to the function Ii in the direction CPj. Then the functions ~i[cpjlsr,lve the system of equations:
This system of equations splits into two independent sets. As mentioned above, the first two equations build up a hyperbolic system. Thus the solvability of the system is guaranteed. The last two equations form a system that is not hyperbolic. The system is converted by adding some regularization terms to the form (5). Then two systems of two equations that can be solved independently are obtained:
4 -10- .6.X3(V)
aX3(V) - Vl- = (17) a Vl a V2 -X4(V) + h(v) -4 aX4(V) aX4(V) - Vl- -10 .6.X4(V) + V2- = (18) a V2 a Vl w
aX3(V) + V2-
2
TVl COSX3(V)
9
+l
sinx3(v)
where the symbol b denotes the Kronecker symbol. The derivative of the functional J with respect to the function Ii in the direction CPj is denoted by J[cpj]' Then the following holds for this derivative:
+ f4(V)
aX5(V) aX5(V) .6.X5(V) + V2- - Vl- = (19) a V2 a Vl -0.IX5(V) + X6(V) + fs(v) -4 aX6(V) aX6(V) -10 .6.X6(V) + V2- - Vl- = (20) a Vl a V2
-10
-4
K
-0. IX 6(V) M
(Vl -
xs)
+ f6(V)
These systems are solved on the domain n = {(x,y) E R2 1x 2 + y2 :S 2 that contains all possible trajectories of the exosystem. Zero Dirichlet boundary conditions are chosen. The border of the domain n is sufficiently far from all the possible trajectories since they are contained in the unit disc. Negligibility of the influence of the boundary conditions on the solution of the system above was verified experimentally.
The implementation of the optimization algorithm will be described next. In order to handle the optimization numerically one has to restrict himself on a much smaller class of the design functions f. The finiteelement-like functions were used. N distinct points (denoted by Pi, i = I, ... , N) in the domain n were chosen. Then continuous functions ; were constructed such that i(Pj) = bi,j. The values of the function i in other points of the set n are determined by a 2D-interpolation. Then the set of the design functions is the set of linear combinations 01 1 + ... + ON N. The design parameters are the coefficients of the linear combinations. Thus the optimum is searched over a finite-dimensional set. Nevertheless the dimension is rather high (2N, since the functions 13 as well as f4 are also to be designed).
The functions 13, f4, fs, f6 on the right-hand side are a kind of design parameters as mentioned above. They are adjusted so that the functional (7) is minimal. This functional splits into two parts that can be minimized independently. The functionals to be minimized are then
The results are demonstrated below. The figure (I) shows the difference of the center manifold values on the trajectory generated by the exosystem Xl = X2, X = -Xl. Solid line represents the values obtained by the finite-element solution of the equation (5) while the dashed line represents the values obtained by the Taylor-series approximation given by Huang in (Huang, 2(00). The oscillations on the solid line are caused by the discretization. They disappear as long as the the discretization is refined. They also influence the behavior of the control loop very weakly. The difference of the finite-element and the Taylor-
Since the minimization problems splits into two parts only the minimization of the functional J 5 ,6 is mentioned. The minimization of the functional J 5 ,6 is carried out by the gradient method.
163
5. CONCLUSION An alternative approach to the nonlinear ou~put·. tracking problem was presented. This approach is based on a direct solution of the partial differential equation using the finite-element method rather than seeking the solution in the form of the Taylor approximations. The advantage of the approach presented here is a possibility of easier algoritmization. Conditions guaranteeing the existence of the solution and convergence of the algorithm are discussed. The approach is demonstrated on the case study of two-cart with an inverted pendulum system.
Fig. 1. The trajectory x(v(t))
REFERENCES C. I. Byrnes, F. Delli Priscoli, A. Isidori and W Kang (1997). Structurally stable output regulation of nonlinear systems. Automatica, 33, 369-285. Devasia, S. (1999). Approximated stable inversion for nonlinear systems with nonhyperbolic internal dynamics. IEEE Trans. onAC, 44,1419-1425. RA. Francis and W Murray Wonham (1976). The internal model principle of control theory. Automatica, 12,457-465. B.A. Francis (1977). The linear multivariable regulator problem, SIAM Journal on Control and Optimization, 15,486-505. J.S.A. Hepburn and WM. Wonham (1981). Error feedback and internal models on differentiable manifolds. IEEE Transactions on Automatic Control, 29, 397-403. Huang, J. (2000). Asymptotic tracking of a nonminimum phase nonlinear system with nonhyperbolic zero dynamics. IEEE Trans. on AC, 45,542-546. Huang, 1. (1995). Output regulation of nonzero systems with nonhyperbolic zero dynamics. IEEE Trans. onAC, 40,1497-1500. Huang J. (2003). On the solvability of the regulator equations for a class of nonlinear systems. IEEE Trans. on AC, 48, 880-885. 1. Huang and WJ. Rugh (1990). On a nonlinear mu1tivariable servomechanism problem. Automatica, 26, 963-972. 1. Huang and WJ. Rugh (1992). Stabilization on zeroerror manifolds and the nonlinear servomechanism problem. IEEE Transactions on Automatic Control, 37, 1009-1013. A. Isidori and C. I. Bymes (1990). Output regulation of nonlinear systems. IEEE Transactions on Automatic Control, 35,131-140. A. Isidori, Nonlinear control systems. Third Edition (Springer, New York, 1995). Roos, H.-G., Stynes, M., Tobiska, L. (1996). Numerical methods for singularly perturbed differential equations. Convection-diffusion and flow problems. Springer, Berlin. Vincent, T.L. and Grantham, W J. (1997). Nonlinear and optimal control systems, New York, WHey.
Fig. 2. The trajectory x(v(t)) - detail
Fig. 3. The difference from the Huang's approach
".:--~-~-----7--~--7---:
Fig. 4. The value of the functional J 5 ,6 series approximations is depicted in the figure (3). The values of the functional (21) are shown in Fig. (4) for five iterations of the minimization procedure.
164
ELSEVIER
PUBLICATIONS www.elsevier.coml1ocnfldifac
UNIVERSAL ALGORITHM FOR THE SOLUTION OF THE REGULATOR EQUATION Branislav Reh:lk .,1 Sergej Celikovsky··,2
• Czech Technical University, Faculty ofElectrical Engineering, Department afControl Engineering, Technickti 2, 166 27 Praha 6, Czech Republic. E-mail: [email protected] .. 1nstitute ofInfonnation Theory and Automation, Academy of Sciences ofthe Czech Republic. P. O. Box 18. 182 08 Prague, Czech Republic. E-mail: [email protected].
Abstract: The computational aspects of the output regulation problem are considered. The standard solution of the output regulation problem uses the explicit solution of the so-called regulator equation being a highly complex PDE combined with algebraic equations. At this moment, the methods to solve regulator equation are, nevertheless, restricted to special classes of systems only. A universal algorithm to solve the regulator equation is proposed here. It is based on a special regularization procedure combin~d ~ith finite elements method and an optimization of the error induced by the regulanzauon. The approach is demonstrated on an example of the two-cart with an inverted pendulum system. Copyright © 2004 IFAC Keywords: Nonlinear control systems, output regulation, partial differential equations.
autonomous exosystem with possible unknown initial states.
I. INTRODUCTION A central problem in control theory and applications is to design a control law to achieve asymptotic tracking with disturbance rejection in a nonlinear system. When a class of reference inputs and disturbances are generated by an autonomous differential equation, this problem is called nonlinear output regulation problem, or alternatively, nonlinear servomechanism problem (lsidori and Byrnes, 1990). The corresponding autonomous differential equation is usually called as the exogenous system. In the sequel, the above setting will be referred to as the output regulation problem. In other words, the output regulation problem treats a possible unknown reference signal and/or disturbances generated by the known neutrally stable
For linear systems the classical output regulation was extensively studied in (Franwon and Murray, 1976; Francis, 1977). For nonlinear systems, the problem was first studied in (Hepburn and Wonham, 1981), and solutions to the output regulation of nonlinear systems have been presented in (Huang and Rugh, 1990; Isidori and Byrnes, 1990) using "full information" which includes the measurements of exogenous signals as well as of the system state. The necessary and sufficient conditions for the existence of a local full information solution of the classical output regulation problem are given in (Isidori and Byrnes, 1990; Huang and Rugh, 1990); they basically mean that the linearized system is stabilizable and there exists a certain invariant manifold. The classical output regulation via error feedback has been solved in (Byrnes et al., 1997; Isidori, 1995) by application of system immersion technique. The plant uncertainty parametrized by
Supported by the Research program J04/98:212300013. Supported by the Gram Agency of the Czech Republic through the research grant No. 10210210709. 1
2
159
unknown constant parameters is treated as a special case of exogenous signals and the solution, extended from the error feedback regulation, is referred to as the structurally stable regulation in (Byrnes et al., 1997).
2. PRELIMINARY FACTS First, ~')me basic facts on nonlinear output regulation problem are recalled (Isidori and Byrnes, 1990; Huang, 2000; Huang, 1995), the notation from (Huang, 2(00) is adopted throughout this paper. Consider the plant
Beginning with (Isidori and Byrnes, 1990), the basic approach to all kind of output regulation problems is to use the explicit solution to the so-called regulator equation which is to be obtained off-line, based on the model of the plant and exosystem only. Nevertheless, regulator equation is a partial differential equation (PDE) combined with algebraic restrictions. Its solvability is a highly complex and questionable issue. From the PDE theory point of view, this is due to its first order structure, not fitting into usual PDE frameworks and, moreover, singularity in its coefficients.
x
= f(x(t)) + g(x(t))u(t) + p(x(t))v(t))
y(t) =
(1)
h(x(t)),
where sufficient smoothness of the vector fields f, g, p and row functions of h is assumed. Further, x(t) E Rn is the state, u(t) E Rm is its input, y(t) E RP its output and w(t) E RV is the so-called exogenous signal. This signal is generated by the so-called exosystem which is supposed to be known and for linear, i.e. for a known (jj x jj )-matrix S and a known (v x jj) matrix Q the exosystem is given by
Existing results solving regulator equation use geometrical approach and special structure of nonlinear systems, for which the solution is trivial. The most simple situation is when overall system (i.e. plant to be controlled plus exosystem) is hyperbolically minimum phase, see already (lsidori and Byrnes, 1990), while (Huang, 2003) shows that for a fairly general class of systems regulator equation may be reduced to its PDE part only. Algorithm to solve this PDE for both minimum and nonrninimum phase cases are suggested in (Huang, 2000; Huang, 1995) based on an undetermined power series technique. Nevertheless, it requires laborious symbolic computations, that are difficult to be implemented as an universal computeraided algorithm applicable automatically to any system from a reasonable class.
iJ = Sv, w(t) = Qv.
(2)
Moreover, the exosystem is assumed to be neutrally stable. This implies the system is Lyapunov stable in both the forward and the backward time direction. Thereby, exogeneous signal is used to describe both reference to be tracked and undesired disturbance to be rejected. This leads to the output regulation problem, which may be tackled by various kind of feedback compensators. The solution to the so-called full information output regulation problem consists in finding the feedback compensator u = a(x, v) such that
The aim of this paper is to provide a universal algorithm for solving regulator equation without referring to some particular structure of the given system. It is based on the regularization procedure that adds new artificial small second order term (Laplacian with small parameter) combined with finite elements methods. The error induced by such a regularization is treated via optimizing certain error functional. The results presented in this paper are encouraging, the proposed algorithm is compared with (Huang, 2000; Huang, 1995) algorithm based on series expansion. It provid es slightly better quality in computing the error zeroing manifold, moreover, this is "automatical" algorithm applicable to any system directly. The theoretical convergence conditions are discussed as well and may be further refined.
(1) the equilibrium x = 0 of the controlled system
is exponentially stable if no exogenous signal is present (2) there exists a neighborhood U C Rn+IJ. of (0, 0) such that for each initial condition (x(O), v(O)) holds limt-++oo e(t) = 0, where the so-called error e(t) is defined as e(t) = h(x(t)) - w(t). The so-called error feedback output regulation problem uses a dynamical feedback compensator using measurements of the error e(t) = h(x(t)) - w(t) only. In both cases, the solution to the following regulator equation is used
8x(v) --a;;-Sv = f(x(v)) + g(x(v))u(v) h(x(v)) = Qv
The paper is organized as follows. Some preliminary facts are presented in the next section, while section 3 presents our algorithm in detail. Section 4 presents illustrative example of the two-cart with an inverted pendulum system. Some concluding remarks are collected in the final section.
(3)
with the condition x(O) = O. More precisely, suppose (x(v), u(v)) is the solution of the regulator equation (3). Then the manifold (x(v), v) is the output zeroing manifold of the augmented plant (1-2) having output h(x) - Qv when with input u = u(v). Based on this property, the feedback compensator solving the full information output regulation problem may be constructed as follows
u = a(x, v) = L(x - x(v))
160
+ u(v),
(4)
bounded domain n c Rn with Lipschitz boundary. The differential equation
L being matrix of gains stabilizing the linear approximation of I. Interested reader may find detailed exposition in (Isidori, 1995), including the error feedback case relying on the solution of regulator equation as well.
8u -a6u + b. 8v has a solution u
3. ALGORITHM FOR SOLVING TIIE REGULATOR EQUATION
V
82
i=1
8v 1
= I: J'i::2
C -
~ div b ~ K > O.
(6)
Corollary: Assume the equation (5) is scalar with k 1 > 0, k 2 > O. Further assume the exosystem (2) is linear and neutrally stable. Then the equation (5) has solution. Proof. The Lyapunov stability means that all eigenvalues have zero real part and none of them is multiple. Hence Trace S = O. Since (using the notation from the lemma above) a = k 1 , b = Sv, C = k 2 we have div b = Trace S = O. Thus the condition (6) is satisfied. D A similar result may be obtained even if the Laplacian 6 is not added. Nevertheless, adding the smallparameter Laplacian is useful since our computational experience shows that the finite-element solution of this augmented equation is "smoother" - exhibits less oscillations on the same mesh without destroying accuracy. Application of the finite-element method for solving the equation is proposed. An extensive treatment of the convergence analysis of numerical schemes for equations like the system (5) is carried out in (Roos, et aI., 1996) and is omitted here. Also, for the sake of brevity investigations of integrability properties of the solution are not presented here. Nonetheless, these properties are crucial for existence of the error functional (see below) as well as for differentiating it. Instead it is assumed all these manipulations are possible.
It consists of two steps - first the control u(v) is chosen and the partial differential equation is solved. Then the control u is adjusted minimizing certain error functional which reflects algebraic part. Then this procedure is repeated until a sufficiently precise approximation of the solution of the system is found. To be more specific, some basic facts on PDE solvability are recalled. Unfortunately, neither existence of its solution nor its uniqueness can be guaranteed for PDE 3. To achieve their solvability their regularization is proposed:
8x(v) -k 1 6x(v) + a:;;-Sv = -k 2x(v) + j(x(v)) + g(x(v))u(v)
.. j, 6
if
a > 0,
The regulator equation (3) has been usually solved using decompositions of the involved functions into the Taylor series in v (Huang, 2000; Huang, 1995). In this paper, such an approach is being replaced by the direct numerical solution of this partial differential equation. As a matter of fact, the system (3) consists of n PDE equations and p algebraic ones in n + m unknown variables x, u. First, without going into details, one may reduce these equation to the case PDE-algebraic system having the same number of equations and variables. The corresponding procedure is similar as the one of (Huang, 2003), nevertheless, not all algebraic equations are to be eliminated in our paper, so that the algorithm of our paper may treat wider class of systems.
+ cu =
(5)
If the error caused by the additional terms in (5) is too large it might be compensated by using the following procedure. Let the symbol r(v) denote the term
with constant vectors k 1 , k 2 with positive (or at least nonnegative) elements.
8x(v) r(v) = a:;;-Sv - (J(x(v))
+ g(x(v))u(v)).
The next issue for (3) is lack of boundary conditions, the only condition posed on the function x is x(O) = O. As it is not possible to solve numerically (5) on the whole space RV, one has to restrict the solution to a sufficiently large bounded domain n c RV and homogenous Dirichlet boundary conditions is used on its boundary. Such a choice of the boundary conditions does not seem to influence the solution significantly if the distance of the trajectory of the exosystem from the boundary of the domain n is large enough and the constants ki are not too large. Thus the condition x(O) = 0 is also satisfied with a fairly good accuracy. Nevertheless, if this condition is unacceptably violated adjusting the boundary conditions is unavoidable.
by adjusting the control u. For the domain holds ~ n. It need not be equal but the set should contain all possible trajectories of the exosystem. The area around the border of the domain n where the boundary condition influence the solution of the differential equation strongly can be out of interest in the error computation.
Theorem 1 «Roos, et aI., 1996) - Lemma 1.18) Let the differential equation in (5) be scalar on a
A formula for the derivatives of the functional J with respect to the control u can be found. The directional
The goal now is to minimize the error functional J =
~
J
(r(v)f r(v)+
(7)
n
(h(x(v) - Qvf(h(x(v)) - Qv)dv
n
161
n n
derivative of the function x(v) in the direction
d dt
X2 X3 X4 Xs
oD x(v)
Sv - kI'~:~.D
x2 h(X2,X3,X.. ,Xs) X4 fs(X2' x3, X4, xs) x6 K -(Xl - xs) M
Xl
~v
x6
(8)
+
(11)
o 1
where Dj, Dg, Dh denote the Jacobi matrices of the functions j, g, h, respectively.
with
h
(9)
can be implemented. Here,
0 is a constant. See (Vincent and Grantham, 1997) for details. In case that the differential equation in (5) does not explicitly depend on the control u this equation can be augmented by adding a right-hand side function as follows:
ox(v)
-k2x(v)+
(2 .
The first two equations represent the dynamics of the cart carrying the pendulum, X I being its position and X2 its velocity. The variables X3 and X4 represent the position and the angular velocity of the pendulum, respectively. The last two variables, Xs and X6 stand for the position and the velocity of the other cart. The variable M = 1.378 means the mass of the cart, m = 0.051 the mass of the pendulum, K = 10 is the spring constant, l = 0.352 is the length of the pendulum, b = 12.98 is the coefficient of friction of the cart, 9 = 9.81 is the gravitational constant, u is the control force.
The simple gradient method in the form
+ a:;-Sv = + \I1(v).
1
. )2 mlx4 smx 3 - bX2 +m (smx3
+bX2 cos X3 - mlx~ sin X3 cos X3 - K(xs - Xl»)
r(vf D
(h(x(v) -
-kl 6.x(v) j(x(v»
= M
-mgcosx3 sinx3 + K(xs - xd) 1 ( . js = l(M + m(sin X3)2) (M + m)g sm X3
The derivative of the functional J in the direction
II
u
o o
oD x(v) D
J
o
- COSX3 l(M + m(sinx3)2)
The fact that these derivatives can be computed will enable to utilize some gradient methods for the optimization. For the derivative of the function r holds:
D
+ m(sin X3)2
M
The zero dynamics is represented by the last four equations. The equations describing the motion of the pendulum build up a hyperbolic zero dynamics while the equations of motion of the added cart are nonhyperbolic zero dynamics.
(10)
We seek a control for the system of equations (11) such that the state Xl tracks the exogenous signal v(t) = Asinwt. This signal is generated by the system
Further steps are analogous as in the previous case. The function \11 is adjusted now instead of the function
u.
iJ
= Sv, S = ( _Ow
~) , Q =
( 1 0) .
The equation of the center manifold (3) attain for our example the form:
4. EXAMPLE: lWO-CART WITH AN INVERTED PENDULUM SYSTEM
V2 0X3 (V) _ V1 0X3 (v) =X4(V) OVl OV2 OX4(V) OX4(V) w2 V2--!:l- - V I - - ! : l - = -l Vl COSX3(V) UVI UV2
The method is demonstrated on the example introduced in (Devasia, 1999) which was extensively studied in (Huang, 2(00). The system consists of two carts which are elastically connected. An inverted pendulum is placed on one cart. The input is the force F acting upon the cart carrying the pendulum. The output is the position of this cart, denoted by Xl. The model of the plant is
TSinX3(V) V2 oxs(v) _ OVl
162
Vl
oxs(v) = X6(V) OV2
(12)
+ (13) (14) (15)
The symbol ~i[cpj] stands for the derivative of the function Xi, i = 5,6 with respect to the function Ii in the direction CPj. Then the functions ~i[cpjlsr,lve the system of equations:
This system of equations splits into two independent sets. As mentioned above, the first two equations build up a hyperbolic system. Thus the solvability of the system is guaranteed. The last two equations form a system that is not hyperbolic. The system is converted by adding some regularization terms to the form (5). Then two systems of two equations that can be solved independently are obtained:
4 -10- .6.X3(V)
aX3(V) - Vl- = (17) a Vl a V2 -X4(V) + h(v) -4 aX4(V) aX4(V) - Vl- -10 .6.X4(V) + V2- = (18) a V2 a Vl w
aX3(V) + V2-
2
TVl COSX3(V)
9
+l
sinx3(v)
where the symbol b denotes the Kronecker symbol. The derivative of the functional J with respect to the function Ii in the direction CPj is denoted by J[cpj]' Then the following holds for this derivative:
+ f4(V)
aX5(V) aX5(V) .6.X5(V) + V2- - Vl- = (19) a V2 a Vl -0.IX5(V) + X6(V) + fs(v) -4 aX6(V) aX6(V) -10 .6.X6(V) + V2- - Vl- = (20) a Vl a V2
-10
-4
K
-0. IX 6(V) M
(Vl -
xs)
+ f6(V)
These systems are solved on the domain n = {(x,y) E R2 1x 2 + y2 :S 2 that contains all possible trajectories of the exosystem. Zero Dirichlet boundary conditions are chosen. The border of the domain n is sufficiently far from all the possible trajectories since they are contained in the unit disc. Negligibility of the influence of the boundary conditions on the solution of the system above was verified experimentally.
The implementation of the optimization algorithm will be described next. In order to handle the optimization numerically one has to restrict himself on a much smaller class of the design functions f. The finiteelement-like functions were used. N distinct points (denoted by Pi, i = I, ... , N) in the domain n were chosen. Then continuous functions ; were constructed such that i(Pj) = bi,j. The values of the function i in other points of the set n are determined by a 2D-interpolation. Then the set of the design functions is the set of linear combinations 01 1 + ... + ON N. The design parameters are the coefficients of the linear combinations. Thus the optimum is searched over a finite-dimensional set. Nevertheless the dimension is rather high (2N, since the functions 13 as well as f4 are also to be designed).
The functions 13, f4, fs, f6 on the right-hand side are a kind of design parameters as mentioned above. They are adjusted so that the functional (7) is minimal. This functional splits into two parts that can be minimized independently. The functionals to be minimized are then
The results are demonstrated below. The figure (I) shows the difference of the center manifold values on the trajectory generated by the exosystem Xl = X2, X = -Xl. Solid line represents the values obtained by the finite-element solution of the equation (5) while the dashed line represents the values obtained by the Taylor-series approximation given by Huang in (Huang, 2(00). The oscillations on the solid line are caused by the discretization. They disappear as long as the the discretization is refined. They also influence the behavior of the control loop very weakly. The difference of the finite-element and the Taylor-
Since the minimization problems splits into two parts only the minimization of the functional J 5 ,6 is mentioned. The minimization of the functional J 5 ,6 is carried out by the gradient method.
163
5. CONCLUSION An alternative approach to the nonlinear ou~put·. tracking problem was presented. This approach is based on a direct solution of the partial differential equation using the finite-element method rather than seeking the solution in the form of the Taylor approximations. The advantage of the approach presented here is a possibility of easier algoritmization. Conditions guaranteeing the existence of the solution and convergence of the algorithm are discussed. The approach is demonstrated on the case study of two-cart with an inverted pendulum system.
Fig. 1. The trajectory x(v(t))
REFERENCES C. I. Byrnes, F. Delli Priscoli, A. Isidori and W Kang (1997). Structurally stable output regulation of nonlinear systems. Automatica, 33, 369-285. Devasia, S. (1999). Approximated stable inversion for nonlinear systems with nonhyperbolic internal dynamics. IEEE Trans. onAC, 44,1419-1425. RA. Francis and W Murray Wonham (1976). The internal model principle of control theory. Automatica, 12,457-465. B.A. Francis (1977). The linear multivariable regulator problem, SIAM Journal on Control and Optimization, 15,486-505. J.S.A. Hepburn and WM. Wonham (1981). Error feedback and internal models on differentiable manifolds. IEEE Transactions on Automatic Control, 29, 397-403. Huang, J. (2000). Asymptotic tracking of a nonminimum phase nonlinear system with nonhyperbolic zero dynamics. IEEE Trans. on AC, 45,542-546. Huang, 1. (1995). Output regulation of nonzero systems with nonhyperbolic zero dynamics. IEEE Trans. onAC, 40,1497-1500. Huang J. (2003). On the solvability of the regulator equations for a class of nonlinear systems. IEEE Trans. on AC, 48, 880-885. 1. Huang and WJ. Rugh (1990). On a nonlinear mu1tivariable servomechanism problem. Automatica, 26, 963-972. 1. Huang and WJ. Rugh (1992). Stabilization on zeroerror manifolds and the nonlinear servomechanism problem. IEEE Transactions on Automatic Control, 37, 1009-1013. A. Isidori and C. I. Bymes (1990). Output regulation of nonlinear systems. IEEE Transactions on Automatic Control, 35,131-140. A. Isidori, Nonlinear control systems. Third Edition (Springer, New York, 1995). Roos, H.-G., Stynes, M., Tobiska, L. (1996). Numerical methods for singularly perturbed differential equations. Convection-diffusion and flow problems. Springer, Berlin. Vincent, T.L. and Grantham, W J. (1997). Nonlinear and optimal control systems, New York, WHey.
Fig. 2. The trajectory x(v(t)) - detail
Fig. 3. The difference from the Huang's approach
".:--~-~-----7--~--7---:
Fig. 4. The value of the functional J 5 ,6 series approximations is depicted in the figure (3). The values of the functional (21) are shown in Fig. (4) for five iterations of the minimization procedure.
164