Integrated Design of Smart Structures Controlled by H∞ Optimal Controllers

European Journal of Control (2003)9:392±406 # 2003 EUCA

Integrated Design of Smart Structures Controlled by H1 Optimal Controllers E. de la Fuente Universidad PoliteÂcnica de Madrid, Escuela Universitaria de IngenierõÂ a TeÂcnia AeronaÂutica, 28040 Madrid, Spain

Control/Structure Integrated Design (CSID) is receiving a significant attention in the field of Smart Structures, specially when related with future Large Space Structures Technology. Robustness of control law in the presence of model uncertainties is a requirement which can be formulated straightforwardly in the frame of H1 control theory. However, its incorporation into a CSID problem is not as easy. The difficulty arises mainly because the H1 controller synthesis must be done iteratively and thus, integration as a part of an optimization scheme requires the development of some additional analytical tools and a careful organization of the optimization routine. This paper presents a general procedure which allows CSID of a Smart Structure linked to an H1 optimal controller.

Control/Structure Integrated Design (CSID) is receiving an increasing attention in the field of Smart Structures, specially when related with future Large Space Structures Technology. Control system designers are used to assume that the plant to be controlled is kept fixed, except perhaps for the placement of sensors and actuators. They design the control law by satisfying a set of requirements expressed as some closed-loop transfer matrix properties. On the other hand, structural designers are not

used to take into account that the structure will possess sensors and actuators and will work in closed loop with a control system, so that its behavior can be modified. Their requirements are placed therefore on the mass, strength, and stiffness of the structure. It is clear that an integrated, i.e. simultaneous design of the structure and the control system, or in other words, simultaneous consideration of control and structural requirements, will necessarily produce a better overall system. This problem is usually posed in the form of a nonlinear optimization problem, with structural and controller parameters taken as design variables and structural and control requirements as objectives. Implementation of control requirements depend heavily on the control synthesis strategy adopted. A variety of them have been treated in the literature. In particular, LQR controllers have been studied in [1±8]. Controllers of LQG type can be seen in [9±11] for instance. Minimax techniques for controller synthesis were used in [12,17]. Controller based in covariance control theory are presented in [18±23]. CSID problems using techniques of independent modal space control are the topic dealt with in [22±26]. Other kinds of controllers based either in Direct Output Feedback, or more or less general Dissipative Controllers have been studied in [27] and [17±30], respectively. H1 optimal controllers, as a part of a CSID problem, have not been studied extensively. In [13], the authors decline to use gradient-based procedures to

Universidad PoliteÂcnica de Madrid, Escuela Universitaria de IngenierõÂ a TeÂcnia AeronaÂutica, 28040 Madrid, Spain. E-mail: [email protected]

Received 6 September 2001; Accepted 3 December 2002; Recommend by S. Weiland and S. Morse.

1. Introduction

393

E. de la Fuente

solve the problem, based mainly on a deficient formulation of the H1 controller existence conditions, stated in terms of the real part of the eigenvalues of the Hamiltonian matrix. In [14], the closed-loop performance and stability are posed in terms of H2 and H1 norms of the closed-loop transfer function, respectively. The optimization problem is performed by genetic algorithms. In [15], the authors deal with a control±structure optimization problem in which constraints are written in the form of H2 and H1 bounds and the structure is assumed to have both parametric as well as unmodelled dynamics uncertainties. An iterative algorithm is proposed in which a selected initial H1-norm bound is taken. The smaller the vallue of for which the controller can be designed, the higher the unmodelled dynamics that can be tolerated. The iterative procedure, that requires the solution of three Riccati equations, is repeated for smaller values of until a satisfactory solution is found, or else the process can be continued until a solution can be obtained for the smallest value of for which a feasible solution can be found. In this case, the minimum upper bound of the H1 norm is found. It will be seen in this paper that a different formulation of the problem allows direct implementation of the H1 optimal controller within the CSID problem. Moreover, the procedure is formulated in a very general way. The interest of the problem dealt with in this paper follows from recognizing that major control requirements are closed-loop performance and stability robustness and that in particular, robustness against model uncertainties is a requirement which can be formulated rather naturally in the framework of H1 control theory, since, for unstructured uncertainties, it can be expressed in the form of constraints on the H1 norm of the closed-loop transfer matrix. However, its incorporation into a CSID problem is not as straightforward. The difficulty arises mainly because the H1 optimal controller synthesis must be done iteratively and thus, integration as a part of an optimization scheme requires the development of some additional analytical tools and a careful organization of the optimization routine. The paper presents a procedure which allows CSID to incorporate requirements stated in terms of H1 norms of closedloop transfer matrix. The paper is structured as follows. Section 2 starts by writing the equations of the open-loop system and formulates the CSID as an H1 optimization problem. The difficulties of solving such a problem are highlighted. Section 3 gives a formulation for the analytical gradient of the infinity norm of a transfer matrix, and the applicability to the problem we are interested

in is given in Section 4. This will lead to the final formulation of the CSID problem as a nonlinear optimization procedure solvable by conventional means. Details of the optimization procedure are given in Section 7. Two examples are finally presented to illustrate the application of the proposed procedure. To ease the reading of the paper, some developments have been left as appendices.

2. Problem Statement 2.1. Mathematical Model of an Active Structure Let us consider a flexible structure equipped with actuators and sensors. Let Mgg ug ‡ Cgg u_ g ‡ Kgg ug ˆ Bgp fp ‡ Bgc fc

…1†

be the dynamic equilibrium equations of the structure. Mgg, Cgg and Kgg are the mass, damping, and stiffness matrices, ug are the structural displacements, Bgp and Bgc are the in¯uence matrices of the perturbation forces, fp and the control forces fc, respectively. Usually, these equations are written in the modal space by means of the transformation1 ug ˆ gh uh in which gh is the modal matrix and uh are the modal displacements, Eq. (1) can be written as a system of differential equations Mhh uh ‡ Chh u_ h ‡ Khh uh ˆ Bhp fp ‡ Bhc fc ,

…2†

where Mhh ˆ Tgh Mgg gh , Chh ˆ Tgh Cgg gh and Khh ˆ Tgh Kgg gh are the modal mass, damping, and stiffness matrices, and Bhp ˆ Tgh Bgp and Bhc ˆ Tgh Bgc are the modal in¯uence force matrices for both perturbation and control forces respectively. For the case of viscous or proportional damping, the above system of equations will be uncoupled, although this condition is not necessary in what follows. Let 0 ug ‡ C1yg u_ g ‡ C2yg ug ‡ vy uy ˆ Cyg

…3† … †

be the sensor output vector, with Cyg the corresponding in¯uence matrices, and vy the sensor noise. The above expression is the most general within the

1 Expression of structural equations in modal coordinates is not necessary in what follows. It is included here because is by far the most common approach in practice.

394

Integrated Design of Smart Structures

limitations imposed by linearity, and may include any type of sensor (displacement, velocity, acceleration, and strain measurement). However, for the sake of simplicity, we will limit ourselves to the case of C2yg ˆ 0, although, as it will be seen in Appendixes B and C, the method presented in this paper is applicable to the general case. Let us finally define a vector ue of structural error variables, i.e. those whose magnitude is intended to be kept as small as possible. Usually these variables will be directly in the form of displacements of some points, or stresses in some places. In any of these two cases, ue can be written in the form ue ˆ Ceg ug :

…4†

Equations (2)±(4) can be written as a ®rst-order linear differential system of equations, in the form

Fig. 1 Block diagram of the smart structure (top) and frequency weighting of inputs and outputs (bottom).

x_ s ˆ Ass xs ‡ Bsw fw ‡ Bsc fc , uz ˆ Czs xs ‡ Dzw fw ‡ Dzc fc , uy ˆ Cys xs ‡ Dyw fw ‡ Dyc fc ,

…5†

where xs ˆ col…uh , u_ h † is the state vector, fw ˆ col (fp, vy) is the total input vector, composed of the perturbation forces fp and the sensor noise. Likewise, uz ˆ col (ue, fc) is a composed error vector, since we will be interested in reducing not only the structural response vector ue,but also the control force vector fc. Otherwise, the problem would not be well posed. Finally, 

 0 I is the state matrix, Mhh1 Khh Mhh1 Chh   0 0 Bsw ˆ is the input in¯uence matrix, Mhh1 Bhp 0   0 Bsc ˆ is the control force in¯uence matrix, Mhh1 Bhc   Ceg gh 0 Czs ˆ is the error vector in¯uence matrix, 0 0     h i 0 0 0 Dzw ˆ , Dzc ˆ , Cys ˆ C0yg gh C1yg gh , Dyw ˆ ‰ 0 I Š 0 0 I Ass ˆ

and Dyc ˆ 0.

 Any input and output may be weighted by appropriate frequency weighting transfer matrices (see Fig. 1). For instance, we can create a new structural response vector uE ˆ WEe(s)ue, in which, WEe(s) is a transfer matrix which is ``large'' only within some frequency range. The effect of taking uE as the new output instead of the original ue is that, in this case, we penalize more severely the structural response ue in that frequency range in which WEe is large. Similar treatment can be applied to any input or output. External perturbation forces may be modelled in frequency giving new vector fW. Finally, dynamics of sensors and actuators may be considered by introducing new weighting frequency transfer functions. Once this augmentation/weighting process is completed, we obtain a system of equations which can be written in the form

2.2. Model Augmentation

x_ S ˆ ASS xS ‡ BSW fW ‡ BSc fc ,

In control theory, the open-loop system equations (5), are usually modified in two respects:

uZ ˆ CZS xS ‡ DZW fW ‡ DZc fc , uy ˆ CyS xS ‡ DyW fW ‡ Dyc fc

 By de®ning additional inputs and outputs (ud, uf), connecting the nominal open-loop system with a number of uncertainty blocks, as is done within the H1 control theory framework, with the intention of representing the ``real'' system. These will increase the length of the input and output error vectors fw and uz,     fw u , uz ) z : fw ) uf ud

…6†

in where uppercase subindexes indicate that the corresponding vector may have been augmented and/or weighted in frequency. 2.3. Remark It shall be emphasized that once the augmenting/ weighting process has been applied (and is to be applied in most practical problems), the specific form

395

E. de la Fuente

of the equations of a structural system (such as those in 5) is destroyed. Thus, the procedure described hereinafter is in fact applicable to any physical system modelled with a linear time invariant mathematical model. 2.4. Controller Implementation To complete the mathematical model of the smart structure, we need to implement the controller. The equations of the controller relate the force control vector fc to the sensor output uy by means of a system of linear differential equations, which in general can be written in the form x_ q ˆ Aqq xq ‡ Bqy uy , fc ˆ Ccq xq ‡ Dcy uy :

…7†

The transfer matrix of the controller is Kcy … j!† ˆ Dcy ‡ Ccq … j!I

Aqq † 1 Bqy :

…8†

The closed-loop system equations are formed by linking (8) together with (6). The ®nal system matrices can be obtained in a straightforward manner, and their expression will not be given for brevity. In the frequency domain, these equations can be written as uZ … j!† ˆ HZW … j!†fW … j!†,

…9†

where HZW is the closed-loop transfer matrix HZW … j!† ˆ DZW ‡ CZS … j!I

ASS † 1 BSW : …10†

Note ®nally that the output and input vectors uZ and fW are composed by (see Fig. 1) 2

3 uE 4 fC 5 uD

2 and

3 fW 4 V 5, uF

respectively. 2.5. Formulation of the Control Structure Integrated Design as an H1 Optimization Problem In a CSID problem, the closed-loop transfer matrix will have two kinds of unknowns. First, a number of parameters p that define the structural design, such as areas, thicknesses, etc. of structural members, and second the matrices that define the equations of the controller, i.e. Aqq, Bqy, Ccq and Dcy, or, equivalently

its transfer matrix, K( j!). In summary, we will have the closed-loop transfer matrix in the form HZW ˆ HZW … j!, p, K†: In H1 control theory framework, we are interested in minimizing the in®nity norm of the closed-loop transfer matrix, i.e. kHZWk1. The augmentation and weighting process described above must be addressed to ensure that the minimization of kHZWk1 will produce an integrated design which is sensible and balanced from the engineering point of view. For the case of unstructured uncertainty, minimization of kHZWk1 will lead to a compromise between stability robustness, measured by the submatrix kHDFk1 and nominal performance, measured by kHEC,Wk1. Thus, the CSID problem can be formulated as a nonlinear optimization problem, in which the objective function is kHZWk1. In addition, other constraints imposed on the structural parameters can be added, such as constraints on mass, stiffness, strength, upper and lower limits of structural parameters, etc. We will join all of these constraints into a constraint vector g(p). In summary, the problem is formulated as min p;K

kHZW …J, , p, K†k1

subject to g…p†  0:

…11†

It is clear that, for a given set of structural parameters (p), for which the structural constraints g(p) are satis®ed, the controller that solves (11) is the H1 optimal controller for the open-loop system corresponding to p. Let us label this controller as Kopt(p). The procedure by which this controller can be obtained will be recalled later on in Appendix B. Thus the problem (11) may be recast in the form min p

kHZW … j!, p, Kopt …p††k1

subject to g…p†  0:

…12†

2.6. Dif®culties Associated with the Solution of the CSID Problem If we attempt to solve (12) by gradient-based optimization techniques, we will find that, since there is no explicit expression for the controller Kopt and that it must be obtained by means of an essentially iterative procedure, the problem would soon become intractable from the economic and numeric point of view. In the following sections, we will present a procedure that successfully overcomes this difficulty and gives a solution to the problem (12). First, we need to develop some mathematical tools. In Section 3, we will

396

Integrated Design of Smart Structures

start by giving a procedure for calculating the gradient of the infinity norm of a general transfer matrix. Afterwards, in Section 4, we will show how to find the gradient of the closed-loop transfer matrix in a way so as to allow the solution to the problem (12). Finally, in Section 5 we will present a final formulation of the problem (12) in a form that is solvable from the mathematical point of view.

3. Gradient of the Infinity Norm of a Transfer Matrix Let x_ ˆ A…p†x ‡ B…p†! z ˆ C…p†x be some generic linear time invariant system. The matrices A, B, and C are assumed to depend on some vector of parameters p. It is assumed that this dependence is known and also that the matrices are differentiable with respect to any element of p. The transfer matrix from w to z will also depend on p, G… j, !, p† ˆ C…p†…J!I

A…p†† 1 B…p†:

Our problem is to ®nd an expression for @ inf …p† , @pk where inf(p) ˆ kG( j!, p)k1 and pk any element of the vector p. Appendix A presents a procedure for calculating the required gradient. It must be recalled that the procedure for obtaining kG( j!, p)k1 gives additionally the frequency !inf at which it is attained. Therefore, both magnitudes are assumed known at the time the gradient is needed. As it can be seen in Appendix A the gradient of the infinity norm of the transfer matrix with respect to a scalar parameter pk can be computed by ! @ inf @G v , …13† ˆ Real uH @p @p k

k !ˆ!inf

where @G=@pk j!ˆ!inf is the gradient of the transfer matrix (which must be of course known) calculated at the frequency ! ˆ !inf, and u and v are the vectors corresponding to the maximum singular value of the transfer matrix G( j!inf, pk). Appendix A shows that, provided the infinity norm is not attained at two or more frequencies

simultaneously, the above procedure will give correct results. Should this not be the case, the infinity norm of the transfer function is not differentiable with respect to pk. This case is analyzed in more detail in Section 6.

4. Gradient of the Closed-Loop Transfer Matrix Coming back to our original problem stated in (12), we find that the expressions developed in the previous section to calculate the gradient of the infinity norm of a transfer matrix, are of no direct application yet, since the transfer matrix referred to in (12) is the closed-loop transfer matrix, and includes the H1 optimal controller, Kopt. To calculate @H/@pk as required by (13), we should need a means to calculate something such as ``@Kopt/@pk.'' However, such calculation is not possible, since Kopt must be obtained through an iterative procedure. To see this clearly, the procedure to obtain the H1 optimal controller is recalled in Appendix B where it is highlighted that the procedure is essentially iterative and therefore cannot be implemented as a part of the optimization problem. The key idea to get a satisfactory solution to (12) is to realize that the Kopt controller is only needed at the end of the optimization process and not necessarily at each intermediate step. In other words, we can approach the optimum values of p and Kopt(p) simultaneously. Only at the end of the process, both shall attain their optimal values. To this end, we augment the vector p of structural design parameters with an additional parameter, which will be called . This new parameter will be used to synthesize a controller, called K, given by Eq. (25), but taking now the scalar  as the value of in (22)± (24). It is clear that the controller K will exist provided that   opt being opt the H1 norm attained with the optimum Kopt controller. The closed loop transfer matrix will now be a function of both p and , i.e. HZW ˆ HZW (p, ). The optimization problem (12) can now be rewritten as min p;

kHZW …s; p; †k1

subject to g…p†  0:

…14†

Remarks. To complete the solution, there are still two remaining dif®culties:  How to calculate the closed loop transfer matrix gradients @HZW/@pk and @HZW/@?

397

E. de la Fuente

This comes from the fact that the controller matrices (25) enter implicitly within the closed loop transfer matrix HZW. However, the gradients can be easily calculated, since now  is a free parameter, and can be treated in the same way as any of p. The details can be seen in Appendix C.  What to do if, in the optimization process, the optimizer requires to calculate kHZW(p, )k1 for a value of  < opt for which there is no solution for the controller K? That is, we must prevent the optimization algorithm to fall into values of  for which there is no solution for the controller. As mentioned, this will happen if the value of  verifies  < opt(p). This precaution is necessary since the optimizer is not informed of that there are forbidden values for . In principle, this would suggest to add a new constraint to (14) of the form  > opt …p† but this would not work, since the evaluation of opt would require a synthesis of the H1 optimal controller and this would lead us to follow the iterative process described above. Instead, since @HZW(s, p, )k1  opt we will write a new constraint in the form  > kHZW …s, p, †k1 …1 ‡ †

…15†

The constraint above could be posed as an equality constraint as well, in which case the optimal Kopt will be certainly obtained (except for the slack, see below) at the end of the process. However, numerical experience acquired by the author shows that the constraint performs properly as an inequality. This is so because the function kHZW(s, p, )k1 is very often a monotonically decreasing function of . Therefore (assuming for the moment that  is zero) the optimizer  are found such will not end until values of p and  that  1 ˆ opt …  ˆ kHZW …s, p, †k p†  which will be precisely the optimum searched. The reason for introducing the slack constant  is that, for the strict optimum, that is, for the closed loop system with the controller Kopt, the closed loop bandwidth is usually very large (even infinite in some cases). From the practical point of view, it is usually interesting to reduce the bandwidth somewhat, and this can be done, at least indirectly, by increasing the value of the parameter  over its minimum value opt.

Besides, the constant  gives an extra protection to the algorithm to prevent incursion into no solution zones for .

5. Final Formulation of the H1 CSID Problem As a result of above considerations, the H1 CSID can be finally written as min p;

kHZW … j!, p, †k1

subject to g…p†  0, …1 ‡ †kHZW … j!, p, †k1

0

…16†

The former (11) or (12) problems that have been shown to be hardly solvable, have now been reformulated as the equivalent (16) problem, which is solvable by standard optimization techniques. It should be stressed out that problem (16) is exactly equivalent to (11) or (12), but solvable by conventional means. This formulation allows for efficient and reliable solution of the H1 CSID problem. However, to see this clearly, the issue of the non differentiability of kHZWk1, as well as the optimization procedure to be followed need separate comments (see the following two sections).

6. Non-differentiability of kHZWk1 Since we are proposing a gradient-based optimization procedure, a natural concern arises about the differentiability of the objective function and the constraints. Although the singular values of a matrix are differentiable in general, this needs not to be the case of the maximum singular value, and that is the case of kHZWk1. This fact is well known and a detailed study can be seen in [16]. Here, we limit ourselves to demonstrate that this issue will not prevent the optimization procedure proposed in this paper to find the right solution. As described in Appendix A, there are two conditions under which kHZWk1 might not be differentiable. These conditions are depicted schematically in Fig. 2. The first one (top) will occur when the infinity norm is attained at two (or more) frequencies simultaneously. This situation could be detected by the optimizer since it will observe abrupt changes in !inf for small changes of the some structural parameter. The second (bottom) will occur when two (or more) singular values of HZW ( j!inf, pk) coalesce into a single one. This situation would lead to significant changes

398

Integrated Design of Smart Structures

w

solution to the quadratic programming subproblem (routine QP). In this case, QP tries to find the best direction of movement, but usually violates the limits imposed to the variables and the structural variables may take negative values, the consequences being catastrophic in our case. The modification consisted simply of moving in the direction suggested by QP but limiting this movement to prevent violation of the limits. Other modifications were needed to avoid endless loops that can be produced occasionally. But the most important modification deals with the strategy to be followed when, despite all the precautions taken, the algorithm falls into forbidden values for . In this case the routine is modified as follows. CONSTR, based on a local quadratic approximation, finds first an optimal direction represented by a vector d, and then calculates the new point xi‡1 from the former one xi by xi‡1 ˆ xi ‡ d

Fig. 2. Two cases in which kHk1 is not differentiable: the infinity norm is achieved at two frequencies simultaneously (top) and coalescence of two singular values of the transfer matrix at !inf (bottom).

in the directions of the unitary vectors associated to the maximum singular value. In any of above cases, one would be tempted to modify the optimization procedure as soon as one of these conditions is detected. This could be done for instance by refining the local approximation, by tracing the hessian eigenvalues, by making local convex approximations, etc. Nevertheless, the numerical experiences performed by the author show that this precaution is in fact unnecessary, and the algorithm, as soon as is ``trapped in the well,'' stops after some few iterations at (or very near to) the optimum point pk.

7. Optimization The algorithm used herein to solve the problem (16) belongs to the sequential quadratic programming category. In fact, subroutine CONSTR from MATLAB [38], designed to solve general non linear constrained optimization problems, has been used with some minor modifications. Some of these modifications were needed to fix some bugs. For instance, the routine usually breaks down when there is no feasible

and the algorithm halves progressively the value of (initially set to one) in an attempt to improve some merit function which is formed as a mixture of the objective function value and the constraint violation. The procedure continues until no further improvement is observed. In the present case, a flag is created to indicate that the value of  is not admissible. This fact is known at the moment of calculating the controller K at the current point, because some of the conditions expressed in Appendix B have failed. If CONSTR finds this flag raised, the merit function value is deliberately modified, forcing the algorithm to halve once more the parameter . The process continues until the value of  moves back into a region of existence of K.

8. Example 1 To illustrate the procedure described above we have chosen a first example, simple enough so as to allow comparison of the results obtained with those of a direct solution. The example consists of three masses m1, m2, and m3 linked by two springs of stiffness k1 and k2. To allow for a meaningful problem, the springs have a mass which is assumed proportional to their stiffness, i.e. mki ˆ ki. We will study two cases (Fig. 3, top)  Case A. The sensor and actuator are both placed at mass 3, and the perturbation force is applied at mass 3. The structural error variable is the displacement of mass 3. In system control jargon, this is a collocated problem.

399

E. de la Fuente Results for Case A 0.82 Sarting point 0.8

| H_ZW | inf

0.78

0.76

Exact 0.74

0.72

0.7 –2 10

–1

0

10

Fig. 3. Mechanical system for the example 1. At the bottom it is shown the multiplicative output uncertainty.

1

10 rho = k2/k1

10

2

10

Fig. 4. Results obtained for Case A. Results for Case B 45

 Case B. The sensor and actuator are placed at mass 2, while the perturbation force and structural error variable are at mass 3. It is therefore a noncollocated problem.

that means that at low frequencies, the model is expected to be accurate within 10%. At high frequencies, the error can be as large as 100%, and at a frequency of 1 rad/s (something about the frequency corresponding to the third mode), the error is 50%. To normalize the outputs ue and fc we will use constant weighting matrices of value 5. This means that values of uE and fC close to unity are acceptable from the engineering point of view, and correspond to ue, fc  5. The complete open loop system has 7 states (6 corresponding to the three normal modes, and the 7th to the uncertainty model), 4 outputs (fC, ud ± the input to the uncertainty block ± and uy) and three inputs ( fW, uf ± the output of the uncertainty block ± and fc). The ®rst three outputs are collected into the vector uZ and the ®rst two inputs into fW. Numerical values are: m1 ˆ 3, m2 ˆ 2, m3 ˆ 1, ˆ 1. The slack variable  is taken as 10 6. Termination tolerance for the objective function and for the design variables were taken as 10 4 and the tolerance for constraint violation was 10 6. Those are precisely the default values of the MATLAB routine CONSTR. The stiffness of the springs are taken as the structural

35 30 | H_ZW | inf

We will assume that the mathematical model is uncertain with a multiplicative output uncertainty (Fig. 3, bottom) represented by the transfer function p s ‡ 2=8 p Wunc …s† ˆ s ‡ 10 2=8

40

25 20 15 1

Exact

10

2

5 0 -2 10

10

-1

0

10 rho = k2/k1

10

1

2

10

Fig. 5. Results obtained for Case B.

parameters, p ˆ [k1 k2]. A constraint is imposed on the mass, so that the total must equal 8 units. In this way, the effect of structural design on closed loop achievable performance is clearly highlighted. This example will show that, by redistributing the stiffness, the actual performance can be significantly improved. The problem to be solved is therefore the following min

k1 ;k2 ;

kHZW … j!, k1 , k2 , †k1

subject to

k1 ‡ k2 ˆ 2,   …1 ‡ 10 6 †kHZW … j!, k1 , k2 , †k1 : …17†

The results are presented in Figs 4 and 5. Although there are two unknowns, k1 and k2, they are however linked by the mass constraint and can be reduced to

400

just one, which is called  and is drawn at the horizontal axis (k1 ˆ 2/(1 ‡ ), k2 ˆ k1). (In the optimization problem however they are considered independent and the mass constraint is explicitly included). In both ®gures, the continuous curve is the exact solution for each value of . They have been drawn by ®rst obtaining the H1 optimal controller corresponding to the plant de®ned by a given value of . To this purpose, the subroutine HINFSYN has been used [38] which follows a procedure similar to that explained in Appendix B. The vertical axis represents the H1 norm of the resulting closed-loop optimal system, and gives the minimum achievable H1 norm for the plant. It is to be noted that for case B, the curve has a pronounced peak (underscaled in the ®gure) at  ˆ 1/3. This is due to the fact that for this value of , the ®rst elastic mode has a node at mass 2, precisely where the actuator and sensor are placed. This means that this mode is uncontrollable and unobservable. The peak observed is therefore simply the open-loop response, which cannot be modi®ed by the controller. In Fig. 4 (case A), the circles represent the progression of the optimization algorithm from an arbitrary starting point. It is seen that the algorithm converges rapidly to the exact solution (marked with an asterisk, although it is hardly visible in the figure). The convergence is attained after 29 iterations. In Fig. 5 (case B), the exact solution shows two minima due to the presence of the peak already mentioned. Starting from point 1, the algorithm converges in 26 iterations to a local minimum and not to the global one. This is completely normal and will happen with all algorithms based on gradient directed search. Starting from point 2, the process converges to the actual global minimum in 51 iterations. In this case, the number of iterations is larger than in the two previous ones, owing to the extremely small curvature at the solution point.

Integrated Design of Smart Structures

Fig. 6. Example 2: Structure and block diagram of the open-loop system.

in closed loop with an H1 optimal controller, and without requiring large control forces. Moreover, this attenuation shall be achieved within some prespecified frequency range. This could simulate for instance the presence of some equipment placed at the middle of the beam and which has some of its natural frequencies within that frequency range. Finally, adequate stability robustness is required against mathematical model uncertainties. The beam has a square cross section. It is assumed to be divided into a number of spans of equal length and constant section. The structural design variables are the sizes of the cross section of each span. The total mass of the beam is constrained to be equal to 10 mass units. The length of the beam is 10 length units. Thus, structural design variables should take values around unity. A viscous damping model is assumed, with a damping coefficient of 0.01, uniform over the frequency. Finally, the Young Modulus is taken as 1000 and the material density is 1000. All data are given in some consistent system of units.

9. Example 2 9.1. Definition of the Problem

9.2. Open-Loop Mathematical Model

The second example illustrates the application of the proposed methodology to a more complex problem and consists of a simply supported beam with variable cross section. The beam is loaded with a perturbation force applied at its middle section and is equipped with two force actuators and two displacement sensors placed at each quarter of the span (see Fig. 6). The overall purpose of the optimization is to keep the displacement of the point of application of the perturbation force as small as possible when working

The open-loop mathematical model is constructed by conventional FEM techniques. Mass and stiffness matrices are built from the corresponding ones of each of the beam elements. Natural frequencies are extracted and the modal viscous damping is added. The dynamic equations of equilibrium are written in the form 2. Finally, the extra inputs and outputs are added and the open loop equations take the form of 5. Of course, these operations are performed at execution time for each structural parameter vector

401

E. de la Fuente

requested by the optimization algorithm. Analytical gradients of natural frequencies and mode shapes with respect to parameters are calculated by conventional means. In particular, the method proposed in [41,42] have been found suitable and have been used in this paper. For the present problem, and according to Fig. 6 the input vector fw consists of the external perturbation force as well as the two outputs coming from the uncertainty model. The output vector uz contains first the displacement of the middle point of the beam, the two control forces, and the two sensor readings uD which will be input to the uncertainty model. Additional inputs to the open loop model are the two control forces fc which will be given by the controller, and additional outputs are the two outputs uF from the uncertainty model, directly injected to the controller. Therefore, the open loop mathematical model has a total of 3 ‡ 2 inputs and 5 ‡ 2 outputs. The lowest three elastic modes are considered in the analysis, so that the number of states is 6.

9.3. Selection of Weighting Factors The problem will now be scaled so as to achieve that kHk1 < 1 corresponds to an acceptable design, being H the closed-loop transfer matrix. To select a reasonable set of weightings, a reference beam of constant cross section of size 1  1 is taken. The static displacement produced by a unit perturbation force is 0.250 length units, and the dynamic displacement is about 12.3. It seems natural to require the controller to eliminate the dynamic amplification, so that a closed loop displacement of the middle point of around 0.250 would be satisfactory. Thus the weighted structural displacement is uE ˆ

1 ue : 0:250

The control forces should be of the order of magnitude of the external perturbation ones, i.e. of the order of 1 force units. Thus fC ˆ fc : The mathematical model of the beam is assumed to be uncertain with a multiplicative output uncertainty such that the error at low frequencies is about 10%. At the frequency corresponding roughly to the third mode (i.e. around 2), the error might be as high as 100%, rising up to 200% at high frequencies. This is

accomplished by a weighting function with a pole at 3.4816 and a zero at 0.17408, i.e. WDd …s† ˆ

8s ‡ 1:3926 : s ‡ 3:4816

The total set of output variables is therefore 2 3 uE 6 7 6 uD1 7 6 7 7 uZ ˆ 6 6 uD2 7 6 7 4 fC1 5 fC2 3 2 32 1=0:250 ue 7 6 76 WYy …s† 6 76 ud1 7 7 6 76 76 ud2 7: ˆ6 WYy …s† 7 6 76 7 6 76 1 4 54 fc1 5 fc2 1 The ``virtual'' equipment placed at the middle point is assumed to have a natural frequency of !n ˆ 0.03 with a viscous damping of 0.1. This has been modelled by assuming that the perturbation force acts on the beam through a second order ®lter of the form p fp …s† !2n ˆ 2 1  2 2 fP …s† s ‡ 2!n s ‡ !2n and values of  ˆ 0.1 and !n ˆ 0.9 have been chosen. Typically, !n would lay around the ®rst natural frequency of the reference beam. The above weighting function produces a maximum force of 1 at a frep quency of !n 1 2 2 ˆ 0:89. To maintain the size of the uncertainty block of the order of 1, the output of the uncertainty block (that will be input to the beam) is weighted as 0.125. Therefore, the total set of inputs to the open loop system is 2 3 fP 6 7 fW ˆ 4 uF1 5 2

uF2

0:1612 6 s2 ‡ 0:18s ‡ 0:81 ˆ6 4

0:125 0:125

32 3 fp 76 7 74 uf1 5: 5 uf2

For the reference beam, the In®nity Norm of the closed loop system, when connected to the H1 optimal controller is found to be 1.17. It is expected that the integrated design will improve substantially this value. The first three elastic modes are considered for the mathematical model of the beam. This makes an

402

Integrated Design of Smart Structures

open-loop system of 10 states (6 for the elastic modes, 2 for the input force frequency weighting and 2 for the uncertainty model). With above weightings, an overall kHk1 of the order of unity is considered to correspond to an acceptable overall design. This would mean that control forces of the order of unity would attenuate the dynamic amplification (i.e. would not produce displacements larger than the static ones) that the beam induces to the equipment in the above specified frequency range, still maintaining acceptable stability robustness margin. As mentioned, the sizes of the cross sections of each span are taken as the structural design variables. The lower and upper limits are taken as 0.1 and 10 length units, respectively. The mathematical formulation of the problem is therefore kHZW …j!, ai , †k1 X subject to a2i li ˆ 10,

min ai ;

  …1 ‡ "†kHk1 , 0:1  ai  10 where ai are the structural design variables,  is the parameter de®ning the controller, " is the slack variable and is taken as 10 4 so that the ®nal controller will be very close to the H1 one. Termination tolerances of the optimization algorithm were taken the same as for Example 1.

specified limits) are carried out. Table 1 shows the results obtained for the final kHZWk1 obtained, along with some statistics. This table shows that, for each case, the problem converges satisfactorily to essentially the same final solution, the standard deviation of the solution being very small. As expected, increasing the number of structural design variables leads to lower values of the final value of the infinity norm of the closed loop system. However, for N ˆ 32, the solution is essentially the same than that obtained for N ˆ 16, so that no further reduction should be expected by increasing the number of spans. The results obtained for the structural design variables can be seen in Fig. 7. The optimal law of variation of the cross section along the span is clearly unsymmetrical. This can be explained easily since a symmetrical variation would have the second elastic mode with a node at the middle section leading to an uncontrollable and unobservable mode. Since the starting point for all the simulations was selected randomly, the problem converges to either cross section laws of the type shown in Fig. 7 or to those mirrored with respect to the middle point. The final value of the Infinity Norm of the closed-loop system is obviously the same in both cases, owing to the symmetry of sensors, actuators and perturbation forces placement. The results indicate that the convergence was correct in all cases, and showed excellent repeatability, even for a large number of structural design variables. The CSID procedure proposed in this paper seemed therefore to work properly.

9.4. Results The problem is solved for a number N ˆ 4, 8, 16, and 32 number of spans, i.e. structural design variables. For each case, ten simulations starting from arbitrary values of the structural design variables (within the

3 32 spans 2.5

16 spans 2 8 spans 4 spans

a(x)

Table 1. Results obtained for the optimum value of kHZWk1 in 10 simulations starting from random initial values of the structural design variables, N is the number of spans taken for each of the simulations.

1.5

1

N 4

8

16

32

kHZWk1 (min)

0.7051

0.6347

0.6156

0.6135

kHZWk1 (mean) kHZWk1 (std. dev.) kHZWk1 (mean)

0.7053 0.0001

0.6348 0.0001

0.6156 0.000

0.6149 0.0021

0.5

0

0

1

2

3

4

5 x

6

7

8

9

10

Fig. 7. Optimal law of variation of the cross section for different number of spans.

E. de la Fuente

10. Conclusions  A gradient-based procedure for performing CSID with H1 control synthesis has been proposed. The procedure overcomes successfully the dif®culties inherent to incorporate H1 controller into the integrated design framework, originated by the iterative procedure which must be followed to synthesize the optimum H1 controller for a given plant, and subsequently the calculation of the gradient.  Two examples are given to illustrate the procedure. For all cases, the procedure proposed in this paper worked satisfactorily.  The procedure is limited to the case of unstructured uncertainties of the open loop mathematical model. Extension to the case of structured uncertainties is possible with some limitations, as well as with some modi®cations of the procedure (see [40]).  Finally, since the method proposed does not make any use of the special form that show the equations of a linear structure, it is also applicable to the integrated design of any linear physical system controlled by the H1 optimal controller.

References 1. Salama M, Garba J, Demsetz L. Simultaneous optimization of controlled structures. Comput Mech 1988; 3: 275±282 2. Slater G, McLaren M. Disturbance model for control/ structure optimization with full state feedback. J Guidance, Control Dyn 1993; 16(3) 3. OÈz H, Khot N. Optimization for efficient structurecontrol systems. Guidance, 1994; 17(6) 4. Khot N, Veley D. Robustness characteristics of optimum structural/control design. Guidance 1992; 15(1) 5. Kajiwara I, Tsujioka J, Nagamatsu A. Simultaneous optimum design of structure and control system to assure stability of higher modes. In AIAA Paper 93-1467-CP, 1993 6. Jin I. Control design variable linking for optimization of structural/control systems. PhD thesis, University of California, 1991 7. McLaren M. Controller methodologies for integrated control/structure design optimization of large space structures. PhD thesis, University of Cincinnaty, 1989. 8. Navid M. Integrated structural and control system design optimization for structures under dynamic loading. PhD thesis, University of California, 1986. 9. Khot N, Heise S. Consideration of plant uncertainties in the optimum structural-control design. AIAA J, 1994; 32(3) 10. Khot N. Integrated structural design and active vibration control. In AIAA Dynamics Specialist Conference, no. AIAA-94±1701-CP, Hilton Head, SC, 1994

403 11. Adamian A. Integrated control/structure design and robustness. PhD thesis, University of California, 1986. 12. Mills R. Robust controller and estimator design using minimax methods. PhD thesis, Standford University, 1992 13. Tsujioka K, Kajiwara I, Nagamutsu A. Integrated Optimum Design of Structure and H1 Control System. AIAA J 1996; 34(1): 159±165 14. Kajiwara I, Nagamatsu A. Integrated Design of Structure and Control System Considering Performance and Stability. In Proceedings of the 1999 IEEE International Conference and Control Applications Hawaii, USA, August 1999 15. Khot NS, OÈz H. Structural-Control Optimization with H2 and H1 norm bounds. Optimal Control Appl Meth 1997; 18: 297±311 16. Giesy DP, Lim KB. H1 norm sensitivity formula with control system design applications. J Guidance, Control Dyn 1993; 16(6) 17. Berg J. Identification and control of lightly damped, large space structures: an experimental evaluation. PhD thesis, Virginia Polithecnic Institute, 1992 18. Skelton R, Ikeda M. Covariance controllers for linear continuous-time systems. Int J Control 1989; 49(5) 19. Hotz A, Skelton R. Covariance control theory. Int J Control 1987; 46(1) 20. Hsieh C, Skelton R, Damra F. Minimum energy controllers with inequality constraints on output variances. Optim Control Appl Meth 1989; 10 21. Yasuda K, Skelton R. Assigning controllability and observability grammians in feedback control. J Guidance 1990; 14(5) 22. Layton J, Peterson L. The benefits of integrated control/structure design using covariance control-a comparison of two approaches. In AIAA 93-1659, 1993; pp. 3126±3136 23. Canfield R, Meirovitch L. Integrated structural design and vibration supression using independent modalspace control. In AIAA Paper 93-1670-CP, 1993. 24. Canfield R, Meirovitch L. Integrated structural design and vibration supression using independent modal space control. AIAA J 1994; 32(10) 25. Shenhar J. Application of control theory to large flexible structures using the Modal-Space control method. PhD thesis, Virginia Polytechnic Institute and State University, 1983 26. Poh S. Modified Independent modal space control of flexible structures. PhD thesis, The Catholic University of America, 1989 27. Lust R, Schmit L. Control-augmented structural synthesis. AIAA J 1988; 26(1) 28. Maghami P, Gupta S, Elliot K, and Joshi S. Integrated controls-structures design methodology: redesign of an evolutionary test structure. J Guidance, Control Dyn 1996; 19: 316±323 29. Fleming P. Computer aided control system design of regulators using a multiobjective optimization approach. In IFAC Control applications of Non linear programming and optimization, 1985; pp. 47±52 30. Padula S, James B, Graves P, Woodward S. Multidiplinary optimization of controlled space structures with global sensitivity equations. Tech. Rep. NASATP 3130, NASA, November 1991

404

Integrated Design of Smart Structures

31. Boyd S, Barrat C. Linear Controller Design. Limits of Performance. Prentice-Hall Inc., 1991 32. Doyle J, Glover K, Khargonekar P, Francis B. State space solutions to standard h-2 and h-inf control problems. IEEE Tran Autom Control 1989; 34(8) 33. Junkins J, Kim Y. Introduction to Dynamics and Control of Flexible Structures. AIAA, 1992 34. Junkins J, Kim Y. First and Second Order Sensitivities of Singular Value Decomposition. J Astronauti Sci 1990; 38(1): 69±86 35. Kwakernaak H, Sivan R. Linear Optimal Control Systems. Wiley-Interscience, 1972. 36. Maciejowski J. Multivariable Feedback Design. Addison-Wesley, 1993 37. Balas GJ, Doyle J, Glover K. Mu Analysis and Synthesis Toolbox for use with MATLAB. The MathWorks Inc., 1993 38. Anon. MATLAB: High Performance Numeric Computation and Visualization Software. Reference Guide. The Mathwork Inc., 1992. 39. Glover K, Doyle JC. State-space formulae for all stabilizing controllers that satisfy an H1 norm bound and relations to risk sensitivity. Systems Control Lett 11: 1988; 167±72 40. de la Fuente E. DisenÄo Integrado de Grandes Estructuras Espaciales. PhD thesis, Polytechnic University of Madrid, 2000 41. Lee, Won, Jung GH. An efficient algebraic method for the computation of natural frequency and mode shape sensitivities. Part I: Distinct natural frequencies. Comput Struct 1997; 62(3): 429±435 42. Lee, Won, Jung GH. An efficient algebraic method for the computation of natural frequency and mode shape sensitivities. Part II: Multiple natural frequencies. Comput Struct 1997; 62(3): 437±443

Appendix A. Calculation of the Gradient of the Infinity Norm of a Transfer Matrix By definition we have [31,36]

inf ˆ max …G… j!, pk ††, !

…18†

where (
) stands for maximum singular value. Since the frequency at which the in®nity norm is attained will be known when the gradient is needed, the above expression may be written

inf ˆ …G… j!inf , pk ††:

…19†

Let us perform a singular value decomposition of G( j!inf, pi) G… j!inf , pk † ˆ UVH , where U and V are complex unitary matrices, and   0 0 ˆ 0 0

with 0 diag (1, 2, . . . , n, being i > 0 the singular values. If we differentiate (19) with respect to p(k), we must take into account that !inf is in general a function of p(k). The gradient of a given singular value of a general complex matrix is (see for instance [33,34])   @i @G ˆ Real uH v i , i @pk @pk where ui and vi are the columns of U and V corresponding to the singular value i. In our case, the above expression reads     @ inf @  @G @!inf @G v : ˆ ˆ Real uH ‡ @pk @pk @!inf @pk @pk Now, the ®rst term     @G @!inf @!inf @G v ˆ v Real uH Real uH @pk @!inf @pk @!inf @!inf @  ˆ 0, ˆ @pk @! !ˆ!inf …20† since, according to (18),  attains a maximum at !inf. Therefore, the expression for the gradient is simply ! @ inf H @G v : ˆ Real u @pk @pk !ˆ!inf

…21†

An obvious equivalent derivation considering the hermitian matrix GHG is ! HG @ inf † T @…G 2 inf ˆ , @pk @pk !ˆ!inf where  is the eigenvector corresponding to the maximum eigenvalue of the hermitian matrix G( j!inf, pk)HG( j!inf, pk)

A.1. Remark From the above derivation, there are two conditions under which the searched derivative might not exist, namely:  The ®nal step in (20) is not correct if @!inf/@pk is not de®ned. This will be the case when, at the point pk, the in®nity norm is attained at two (or more) frequencies !inf 1 and !inf 2 simultaneously. In this case, the function !inf (pk) is not continuous.

405

E. de la Fuente

 The frequency !inf is continuous at pk, but the maximum singular value of the G is not single at pk. In this case, the unit vectors u and v will not be continuous at pk. Under any of the above conditions, the existence of @kGk1/@pk cannot be assumed.

where and X and Y are the solution of the Algebraic Riccati Equations ATSS X ‡ XASS X…BSc R 1 BTSc ASS Y ‡ YATSS

Y…CyS V 1 CyS

2 W†X ‡ Q ˆ 0,

2 Q†Y ‡ W ˆ 0

and the matrix Z is given by

Appendix B. Synthesis of H1 Optimal Controller We briefly summarize the procedure to synthesize the optimal H1 controller for a given plant. Further details can be seen in [31,32,36]. We start from an open loop plant described by the system of state space Eq. (6). First, we will assume that the open loop plant verifies the following relationships: DTZc CZS ˆ 0, BSW DTyW ˆ 0. It is usual in practice to assume no cross-weighting of error variables and control forces, as well as no cross-coupling between perturbation forces and sensor noise. However, as shown in ([36]), the general system (6) can always be cast in this form, by simple input and output scaling. We will assume that this has already been done. Besides, for the sake of simplicity we will further assume that DZW ˆ 0 and Dyc ˆ 0 (see the comment at the end of this appendix). Finally, we will assume that     ASS J!I BSW ASS J!I BSc and CyS DyW CZS DZc are respectively of full column and full row rank for all real Z. Concerning the two last conditions, the ®rst one implies that there are not frequencies at which the actuators do not excite the Z variables, and the second means that there is not any frequency at which the external perturbations do not excite the sensors. It is easy to verify that both conditions are ful®lled automatically in our present formulation since the variables Z contain indeed the actuator forces fc and the perturbations do include the sensor outputs. Furthermore, we will call R ˆ DTZc DZc , T T T V ˆ DyW DyW , W ˆ BSW BSW and Q ˆ CZS CZS , where R and V must be positive definite for the problem to be well posed. Now, following Doyle [32], given a real positive number , there exists a controller K such that kHZWk1 < if and only if  X ˆ XT  0,  Y ˆ YT  0,  Z ˆ ZT  0,

Z ˆ X…I

2 YX† 1 :

…24†

Furthermore, if these conditions are satis®ed, one controller (in fact, there are in®nite) that achieves kHZWk1 < is given by Ccq ˆ R 1 BTSc Z, Bqy ˆ YCTyS V 1 , Aqq ˆ ASS

BSc Ccq

Bqy CyS ‡ 2 YQ:

…25†

This speci®c controller is often called the minimum entropy controller [31], since it minimizes the so called

-entropy of the transfer matrix HZW. The procedure to find the controller K1 that minimizes kHZWk1 is clearly iterative, and is performed by a bisection algorithm. It is now clear that, as mentioned in the paper, the above controller synthesis is iterative and so could hardly be implemented as a part of the optimization problem (12). For the general case, that is DTZC CZs 6ˆ 0, BSW DTyW 6ˆ 0, DZW 6ˆ 0, the expressions of the controller matrices become considerably longer, but in any case, for a given value of , they can be obtained explicitly (see for instance [36] or [39] for the details). This will be the case, for instance, if accelerations are taken as sensors and/or as output error variables.

C. Gradient of the Closed-loop Transfer Matrix We will assume that the open-loop state matrices depend on a vector of parameters p. Furthermore, the controller is parameterized with an extra parameter . Therefore, the closed-loop transfer matrices will be dependent on both p and . Once this dependence is explicitly known, the gradient of the Infinity Norm of the closed-loop transfer matrix can be calculated by using (13). The purpose of this appendix is to illustrate the calculation of the gradient of state space controller matrices with respect to both pk and , where pk can be any element of p. From these results, the gradient of the closed-loop matrices will be immediate. The

406

Integrated Design of Smart Structures

following lines illustrate the procedure to follow for the simplified (although very common) case detailed in Appendix B. At the time the gradient is searched for, the solution X and Y of the two Riccati equations, as well as the infinity norm of the closed-loop transfer matrix, inf, are known for the current values of p and . Also, it is assumed that the gradient of the open-loop state matrices with respect to pk are known. With these preliminaries, the calculation of the gradient of the controller matrices proceeds as follows. Differentiating both (22) and (23) we get two linear Lyapunov equations in X and Y, namely AT1 X ‡ XA1 ‡ Q1 ˆ 0,

…26†

A2 Y ‡ YAT2 ‡ W2 ˆ 0

Once both X and Y are known by solving the Lyapunov equations (26) we can obtain the variation of the controller matrices very easily. First,

 … 2 3 YX

A1 ˆ ASS

…BSc R

A2 ˆ ASS

…CyS V

2

Ccq ˆ

‡ R 1 BTSc Z, Bqy ˆ YCTyS V

BSc Ccq

1

BSc Ccq

 2 Q†Y,

 2 YQ

 2 YQ:

 2 W†X ‡ Q,  2 Q†Y ‡ W: 1

CyS

 2 W† ˆ BSc R 1 BTSc ‡ BSc R 1 BTSc BSc R 1 RR 1 BTSc 2 3 W

…CyS V 1 CyS

‡ YCTyS V

Bqy CyS

The terms …BSc R 1 BTSc  2 W) and (CySV  2Q) can be obtained from …BSc R 1 BTSc

1

Bqy CyS

W2 ˆ ASS Y ‡ YATSS Y…CyS V 1 CyS

 2 YX†

R 1 RR 1 BTSc Z ‡ R 1 BTSc Z

 W†X,

Q1 ˆ ATSS X ‡ XASS X…BSc R 1 BTSc

 2 YX

1

and now

Aqq ˆ ASS BTSc 1 CyS

 2 YX†

1

 2 YX†

 …I

…I

YCTyS V 1 VV 1 ,

with 1

1

 2 YX†

Z ˆ X…I

 2 W,

 2 Q† ˆ CyS V 1 CyS ‡ CyS V 1 CyS CyS V 1 VV 1 CyS 2 3 Q

 2 Q:

2 3 YQ

and we get ®nally in explicit form the variation of the controller matrices in terms of the variations p of the open-loop plant matrices with p, and the variation  of the extra parameter . The variation of the system closed-loop matrices is now immediate. For the general case, that is DTZC CZs 6ˆ 0, BSW DTyW 6ˆ 0, DZW 6ˆ 0, DyC 6ˆ 0, the procedure follows the same steps as before since first the gradient open-loop matrices with respect to pk can be calculated exactly in the same manner, and second, the controller matrices are explicitly known in terms of the parameter  as mentioned in Appendix B. The resulting expressions are considerably more involved, but there is no inherent difficulty to arrive at them.