Gain scheduling based on linear parameter-varying (LPV) systems has been a popular method for nonlinear control system design in the past decade^[1]. Using linearization techniques,the inherently nonlinear system is described as a set of linear models,which are driven by several so-called scheduling parameters that can be measured online or a-priori known. Then based on these LPV models,the controller is synthesized in a linear-like fashion.

For LPV systems with uncertainties,controller synthesis to achieve optimal robustness is a challenging problem,as both scheduling parameters and uncertainties can enter the model in arbitrary ways. Several robust techniques were applied to uncertain LPV problems, ${\rm e.g.,}$ $\mu $-synthesis^{[2, 3]} and Lyapunov-based methods^{[4, 5, 6]},all of which involved solving some linear matrix inequality (LMI) based convex programs. The parameter dependence of the LMIs is not easy to handle,one natural way was to grid the parameter set^[4],which,although is a direct way,cannot guarantee the feasibility of solution on the entire parameter set. One way bypassing gridding was to restrict the dependences on either uncertainties or scheduling parameters to some specific classes, e.g.,affine structure^[5],and linear fractional transformations^{[2, 3, 6]},making problems more tractable while bringing about some conservativeness.

Probabilistic robust control (PRC) alleviates the drawbacks of the methods mentioned above by shifting the usual deterministic sense of robustness to a probabilistic one. A controller that is robust for most of uncertainty samples is satisfactory,which means a risk of instability is acceptable if it is very small. This evolutionary concept makes many robust design problems easier to compute. Moreover,PRC methods can handle very general uncertainty structures compared to the existing deterministic approaches,and the conservativeness of results is controllable.

In fact,PRC has received lots of attentions in control community since late 1990s. Various algorithms and their applications are documented,which are classified into three main methodologies by Calafiore et al.^[7]: the approach based on the statistical learning theory,the sequential approximation methods based on gradient iterations or stochastic ellipsoid (SE) iterations,and the scenario approach (SA). SE algorithms^{[8, 9, 10]} proceed in an iterative scheme,online generating scenarios for oracle,and the sample size per iteration is independent of both the number of uncertainty parameters and the number of the design variables. So SE approaches seem to perform quite efficiently,and have been the most popular PRC methods currently. However,for an SE algorithm,the total amount of scenarios consumed by all the iterations is very large,even much larger than that of SA. So for controller synthesis problems of LPV systems,whose linearization procedure may need a great deal of calculation,SE methods might be too time-consuming.

SA,developed by Calafiore and Campi^{[11, 12]},gives a one-shot solution to a convex optimization problem with constraints generated through uncertainty randomization. The number of samples that ought to be made in the randomization procedure is a first-order polynomial to the number of design variables^[13]. So when facing high-dimensional design problems,such as robust semidefinite programming,the computational load of SA is very heavy,and even the computing resources would be used up. This is why SA is more valuable in theory than in practice.

The main contribution of this paper is three-fold. Firstly,we propose a batch iterative algorithm for SA,which demonstrates a good solution performance as the SE methods for LPV controller design,while being more efficient. Secondly,in the probabilistic robust LPV control context,using randomization technique,we provide a unified and concise way to cope with the parameter dependency of the LMI constraints,thus remarkably reducing the number of constraints needed for optimization computation,compared to the commonly-used gridding approach^[1]. Thirdly,the proposed method is applied to design a robust LPV controller for a small helicopter. To the best of our knowledge,small helicopters would be one of the most difficult aerial vehicles to control. The state of art has explored various control synthesis methods,focusing on the nonlinear and robust design issues (see,for example,Wang et al.^[14],Liu et al.^[15],Cai et al.^[16],and Simplicio et al.^[17]). Both problems of robustness and nonlinearity are dealt with simultaneously in this paper,and such a challenging task is a true touchstone to the proposed controller synthesis algorithm.

The paper is organized as follows. In Section II,the scenario design problem is introduced,where two iterative algorithms, including the batch one,are proposed. Section III gives a sufficient condition for the existence of a solution for a robust output-feedback LPV control problem,where the controller formulation is presented,and computational issues associated with iterative scenario approach (ISA) are discussed. In Section IV, using the batch algorithm developed in previous sections,a robust output-feedback controller is designed for a small helicopter with uncertainties in multiple physical parameters,and performance is verified through both comparisons and nonlinear simulations. Conclusion is made in Section V.

A general robust convex program (RCP) takes the following form:

where ${{\Delta }} \subseteq {{\bf R}^{{n_{\pmb{\delta }}}}}$ is the uncertainty set and $\varphi ({\pmb{\theta }},{\pmb{\delta }} )$ is convex in ${\pmb{\theta }} \in {{\Theta }} \subset {{\bf R}^{{n_{\pmb{\theta }} }}}$. Despite the convexity assumption on both objective and constraints,RCP is in general an NP-hard problem due to the possible infinity of constraints.

SA looses the constraint satisfaction from ``all" to ``almost" by applying the so-called randomization technique on the uncertainty parameters. Assume that $N$ independent identically distributed samples of ${\pmb{\delta }},\;\left\{ {{{\pmb{\delta }} ^{(\rm 1)}},{{\pmb{\delta }} ^{(\rm 2)}},\cdots,{{\pmb{\delta }} ^{( N)}}} \right\}$,named ``scenarios",are extracted according to some probability $\Pr ({\pmb{\delta }} \in {{\Delta }} )$,then a scenario convex program (SCP) is constructed as

Solution to problem (2) guarantees the robustness in a probabilistic sense,which is stated by the following theorem.

Theorem 1^[13]. Fix two real numbers $\varepsilon \in (0,1)$ (level parameter) and $\beta \in (0,1)$ (confidence level). If

then with probability of no smaller than $1 - \beta $,one of the following two things happens: 1) ${\rm{SCP}}_{{N}}$ is unfeasible, and then (1) is unfeasible; 2) ${\rm{SCP}}_N$ is feasible and then its optimal solution ${{{\hat {\pmb{\theta }} }}_{N}}$ is $\varepsilon $-level robustly feasible,which means $\Pr \{ {\pmb{\delta }} \in {{\Delta }} :\varphi ({\pmb{\theta }},{\pmb{\delta }} ) > 0\} \le \varepsilon $. Although a lower bound of the constraint number has already been given,the computational load of a one-shot solution to ${\rm{SCP}}_N$ is still too heavy. For example,if we set $\varepsilon = 0.05,\;\beta = {10^{ - 5}},\;{n_{\pmb{\theta }}} = 157$,as seen later in Section IV,then 6742 constraints are generated according to (3). This solving scale exceeds the capability of most current convex optimization softwares.

Inspired by Wada$'$s work^[18],we present here an iterative scheme based on the maximum volume ellipsoid (MVE) cutting plane method,which is known to be particularly efficient for large scale convex optimization problems.

Firstly,we introduce the MVE within a polytope. It is assumed that the parameter ${\pmb{\theta }}$$'$s space is an $ n_{\pmb{\theta }}$-dimensional Euclidean space. Then,an ellipsoid in the parameter space is defined as

where ${Q} = {{Q}^{\rm{T}}} \in {\bf{R}}_{\rm{ + }}^{{n_{\pmb{\theta }}} \times {n_{\pmb{\theta }}}}$ defines the ellipsoid$'$s shape, ${{\pmb{\theta }}_{{c}}}$ is the ellipsoid center. The volume of the ellipsoid is given by ${\mathop{\rm vol}\nolimits} {\rm{(}}E({ Q},{{\pmb{\theta }}_{{c}}}{\rm{)) = det(}}{Q}{\rm{)}}$. Given a convex polytope,${{P}}= \{ {\pmb{\theta }} \in {{\bf{R}}^{n}}:{\pmb a}_i^{\rm{T}}{\pmb{\theta }} \le {{b}_i},i = {\rm 1},\cdots,L\} $ the MVE within a polytope $P$ is defined by

It is easily seen that the MVE can be found by solving the following optimization problem:

To facilitate the plane-cutting procedure,we limit the search space to a hypercube that is known a priori.

Assumption 1. Let ${{\pmb{\theta }}_1}$ be the initial guess of parameter ${\pmb{\theta }}$,and define ${{\pmb{r }}_m} = {\left[{{{ r}_1}\;{{ r}_2} \cdots {{r}_{{{n}_{\pmb{\theta }} }}}} \right]^{\rm{T}}}$,then the hypercube ${{{\Theta }} _1} \subset {{{\bf R}}^{{n_{\pmb{\theta }} }}}$,which is assumed to contain the optima ${{\pmb{\theta }}^*}$of problem (11),is given as

where ${{r}_i} > {\rm{0}}$,and ${{{\pmb e}}_i} \in {{\bf{R}}^n}$ is the unit vector with its $i$th element being 1.

Then a basic iterative algorithm to solve problem (2) is presented as follows.

Algorithm 1.

Inputs: ${{\pmb{\theta }}_1}$; ${{\pmb r}_{\rm{m}}}$; $\{ {{\pmb{\delta }} ^{(k)}}\} _{k = 1}^N$; $\tau $.

Initialize: $l = 1$,${{{A}}_1} = {\left[{{{\pmb{I}}_{\rm n}}\;\; - {{\pmb{I}}_{\rm n}}} \right]^{\rm{T}}}$,

${{{B}}_1} = {[{{\pmb{\theta }}_1} + {{\pmb{r}}_{\rm{m}}}\quad - {{\pmb{\theta }}_1} + {{\pmb{r}}_{\rm{m}}}]^{\rm{T}}}$,$\{ {{\bf{\hat {\pmb{\theta }} }}_j}\} = {\rm{NULL,}}~~ j = 0$ where $j$ counts the feasible solutions. The initial polytope,which is a hypercube,is defined as

1) Perform oracle on the $N$ scenarios.

If $\varphi ({{\pmb{\theta }}_l},{{\pmb{\delta }} ^{(k)}}) \le 0$ for all $k \in \left\{ {\rm 1,2},\cdots,N \right\}$ ,then a feasible solution is obtained,and we set

else set

2) Update. Add a new cutting plane to ${P_l}$ :

find the MVE ${E_{l + 1}}$ in the new polytope:

and obtain $({{Q}_{l + 1}},{{\pmb{\theta }}_{{c}}})$. Set ${{\pmb{\theta }}_{l + 1}} = {{\pmb{\theta }}_{{c}}}$.

3) Stopping rule. If ${\left\| {{{\pmb{c}}^{\rm{T}}}\left( {{{{\bf{\hat {\pmb{\theta }} }}}_j} - {{{\bf{\hat {\pmb{\theta }} }}}_{j - 1}}} \right)} \right\|_2} \le \tau $,the algorithm stops, else set $l = l + 1$,and go to step 1).

Outputs: ${\hat {\pmb{\theta }} _j}$}.

The algorithm presented above consists of two iteration loops. The inner is for feasibility verification,while the outer is for update. Notice that the number of linear inequalities that comprise ${{P}_l}$ increases by one for each outer loop iteration. Consequently,the growing computational burden would cause numerical instability,especially for some high dimensional cases. Therefore a mechanism to remove the redundant ones from the constraint set is quite necessary. In this paper,a heuristic constraint-pruning logic,suggested by ^[19],is adopted. We fix the number of linear inequalities in ${{P}_l}$ to be an integer between $3{n_{\pmb{\theta }}}$ and $5{n_{\pmb{\theta }}}$ typically,and then rank the constraints according to the relevance ratio,defined by

where ${{F}} = {n_{\pmb{\theta }}}{{Q}_{l + 1}}$ defines the shape of an ellipsoid centered at ${{\pmb{\theta }}_{{c}}}$,which is supposed to cover the polyhedron ${{P}_{l + 1}}$. And then we substitute the least relevant inequality (the inequality with the largest relevance ratio) with the new one. Thus (5) of the update procedure in Algorithm 1 is modified as

where ${{\pmb{a}}_n}$ is the $n$th row vector of ${{{A}}_{l + 1}}$,and ${b_n}$ is the $n$th element of ${{{B}}_{l + 1}}$. In this paper,we limit the total number of inequalities to be $3{n_{\pmb{\theta }}}$.

The convergence of the algorithm is guaranteed by the characteristic of the MVE cutting plane method. A theoretical upper bound for the number of outer iterations of the MVE-based SE method,suggested in ^[18],is also applicable to Algorithm 1,despite its conservativeness. For details about cutting plane methods,readers can refer to ^[19]. If Assumption 1 is not satisfied,there would be otherwise two special cases. First,if the primal ${\rm{SCP}}_N$ problem (2) is unfeasible,the algorithm will not output any feasible results after proceeding for a large quantity of cycles, even if the number of iterations reaches its upper bound. Second,if ${{\Theta_1 }}$ does not contain the optima ${{\pmb{\theta }}^*}$, which is possibly encountered in practice as the initial hypercube may fail to satisfy Assumption 1,the algorithm just finds a sub-optimal solution within ${{{\Theta }}_1}$.

Given an ${\rm{SCP}}_N$ problem,Algorithm 1 proceeds in a deterministic framework,although the underlying interpretation is in a probabilistic sense. Different from SE methods,the scenarios used for test are prepared in advance. However,for Algorithm 1, it is still cumbersome to test all the $N$ scenarios in every inner iteration,because the quantity $N$ depends on the dimension of the optimization parameter (see (3)),usually being a large number,larger than that used in SE methods. Consequently, Algorithm 1 has little advantage over the SE methods in respect of algorithm$'$s total execution time.

The basic idea to amend Algorithm 1 is to reduce the number of scenarios consumed in each inner loop. One natural way is to randomly batch the $N$ scenarios for verification; the underlying idea is that successive demonstrations of feasibility for the randomly-sampled batches will steadily enhance our confidence on the solution$'$s feasibility for all the $N$ scenarios. To illustrate the rationale of this idea,we present the following theorem.

Theorem 2. Suppose we have to inspect a set of $N$ scenarios. A $k$-length subset randomly sampled from the $N$ scenarios,called a $k$-sample batch,is defined to be ``good'' if and only if the solution is verified to be feasible for all the $k$ scenarios. Then,if $m$ successive $k$-sample batches have been verified to be good,our confidence that the solution is feasible for all the $N$ scenarios,denoted by

Proof. For a given solution,let $\lambda $ be the chance of a scenario to be good,then in a $k$-sample batch,the number of good scenarios,denoted by ``$x$'',obeys the binomial distribution,i.e.,

where $l\left( {x|\lambda } \right)$ is a likelihood function. If $x = k$,(9) is simplified as $l\left( {k|\lambda } \right) = {\lambda ^k}$.

Then for $m$ successive good $k$-sample batches that are $\rm i.i.d$,the likelihood function can be written as

In the Bayesian sense,the posterior probability of $\lambda $ is given by

where $g(\lambda )$ is $\lambda $$'$s prior distribution. For a given $N$-length set of scenarios,$\lambda $ is a discrete stochastic variable that can be assigned $(N + 1)$ values: $0,1/N,2/N,\cdots,\;1$. Once a good $k$-sample batch has been obtained,it is definitely known that $\lambda \ge k/N$. And further,from a Bayesian point of view,if we have no preference about $\lambda $$'$s evaluation in interval $\left[ k/N,\;1\right]$,it is most reasonable to assume $\lambda $ to be uniformly distributed in this interval. So $g(\lambda )$ is given as

Inserting (10) and (12) into (11),the confidence that $\lambda $ equals to 1 is then obtained as

Theorem 2 indicates the number of times we have to repeat for sampling inspections if $k$ and $\eta $ are determined in advance. It then inspires Algorithm 2,for which we just give the iteration part below,as the rest of Algorithm 2 is similar to Algorithm 1.

Algorithm 2.

Iterations:

1) Randomly pick $K$ samples from the set of $N$ scenarios,and perform oracle on this subset.

If $\varphi ({{\pmb{\theta }}_l},{{\pmb{\delta }} ^{(k)}}) \le 0$ for all $k \in \left\{ {1,2,\cdots,K} \right\}$

If ${\rm{Do\_Batch\_Check }} \ne 1$,set

else set $fcnt = fcnt + 1,\quad {{\pmb{\theta }}_l} = {{\pmb{\theta }}_{l - 1}}$,and go to step 3);

else set $f = \varphi ({{\pmb{\theta }}_l},{{\pmb{\theta }} ^{(k)}}),\;{\pmb{g}} = {\partial _{\pmb{\theta }}}\varphi ({{\pmb{\theta }}_l},{{\pmb{\delta }} ^{(k)}}),\;fcnt = 0.$

2) Update. Add a new cutting plane to ${P_l}$ according to (7),find the MVE ${E_{l + 1}}$ in polytope:

obtain $({{Q}_{l + 1}},{{\pmb{\theta }}_{{c}}}),\;$ and set ${{\pmb{\theta }}_{l + 1}} = {{\pmb{\theta }}_{{c}}}$.

If ${\left\| {{{\pmb{c}}^{\rm{T}}}\left( {{{{\bf{\hat {\pmb{\theta }} }}}_j} - {{{\bf{\hat {\pmb{\theta }} }}}_{j - 1}}} \right)} \right\|_2} \le \tau $,set ${\rm{Do\_Batch\_Check = 1}}$.

3) Stopping rule.

If $fcnt > M$,the algorithm stops,

else set $l = l + 1$,and go to step 1).

In the algorithm,$N,K$ and $M$ respectively stand for the total amount of scenarios generated by the SCP problem,the number of scenarios comprising one batch and the number of successive good batches.

Remark 1. For Algorithm 2,at the beginning,everything goes the same way as in Algorithm 1 except the number of scenarios used for feasibility verification,until the objective function ${{\pmb{c}}^{\rm{T}}}{\pmb{\theta }}$$'$s descending rate reaches a prefixed threshold $\tau $. Then the algorithm stops improving ${{{\pmb c}}^{\rm{T}}}{\pmb{\theta }}$,and tries to find a solution that is able to pass successive verifications on $M$ randomly-picked $K$-sample batches.

Remark 2. Only the last one of the output set $\{ {{\bf{\hat {\pmb{\theta }} }}_j}\} $ is the solution we want,and we recommend a complete feasibility verification of this solution on all the $N$ scenarios. Fortunately,if we have set a high confidence,e.g.,$\eta = 0.95$,the last one is almost always the ``right" one.

Remark 3. The above batch algorithm is designed for robust optimization problems,and a slight modification to the algorithm will fit it for robust feasibility problems,while merely removing the parts about objective function improvement away from the iteration process.

Given two compact sets ${\cal P} \subset {{\bf{R}}^{{n_\rho }}}$ and ${{\Delta }} \subset {{\bf{R}}^{{n_{\pmb{\delta }} }}}$, consider an uncertain LPV system ${{\bf{\Sigma }}_{{\cal P},{{\Delta }} }}$,i.e.,

where ${\pmb{x}} \in {{\bf{R}}^{{n}}},{\pmb{z}} \in {{\bf{R}}^{{{{n}}_{\pmb z}}}},{\pmb{y}} \in {{\bf{R}}^{{{{n}}_{{\pmb y}}}}},{\pmb{d}} \in {{\bf{R}}^{{{{n}}_{{\pmb d}}}}},{\pmb{u}} \in {{\bf{R}}^{{{{n}}_{{\pmb u}}}}}$ represent the states,performance outputs,measured outputs,disturbances and control inputs, respectively; $\rho \in {\cal P}$ is the scheduling parameters, which is a-priori unknown but can be measured online; and ${{\delta }} \in {{\Delta }} $ denotes the uncertainties in model parameters. For succinct expression,we omit the dependence on time of the variables and parameters in (13). It is known that ${{\pmb{\Sigma }}_{{\cal P},{{\Delta }} }}$ can be derived from the Jacobi linearization of a nonlinear model. An output feedback LPV controller ${K_{\cal P}}$ scheduled by $\rho $,can be written as

In this paper,we are concerned with the induced L2-norm performance of closed-loop LPV system (13) and (14),which is defined as

where ${F_{\cal P}}$ denotes the set of all piecewise continuous functions mapping time $t$ into ${\cal P}$ with a finite number of discontinuities in any interval^[20]. Following the loop-shifting techniques^[21],it is possible to transform any general LPV plant as in (13) to a simplified form $\sum\nolimits_{{\cal P},{{\Delta }} }^S {} $,written as (15), through norm-preserving transformations on ${\pmb z}/{\pmb d}$ and invertible transformations on ${\pmb y}/{\pmb u}$:

For the L2-norm robust control of an LPV system in the presence of model uncertainty,we give a sufficient condition as follows

$\sum\nolimits_{{\cal P},{{\Delta }} }^S {} $,and a controller ${K_{\cal P}}$,the closed-loop system is guaranteed to have an L2-norm performance level $\gamma > 0$ if there exist two symmetric matrices ${{X}} > 0$ and ${{Y}} > 0$,such that for all $\rho \in {\cal P}$ and ${\pmb{\delta }} \in {{\Delta }} $,the following LMIs hold:

where

and

The theorem can be easily proved by directly applying Wu$'$s theorem (see ^[20],Theorem 3.3.1) on parameter dependent output-feedback controller synthesis,so long as we treat uncertain parameter ${\pmb{\delta }} $ the same way as scheduling parameter $\rho $,and use the single quadratic Lyapunov function (SQLF) as a degraded form of the parameter-dependent Lyapunov function (PDLF). Once such ${{X}}$ and ${{Y}}$ are found,we can then construct the robust output-feedback LPV controller based on the nominal state-space matrices. Define^[20]

where

and

And further,let

then the $n$-order strictly proper controller ${K_{\cal P}}$ achieving $\gamma$-performance is presented as

The RCP problem associated with Theorem 3 can be expressed as

where ${\pmb{\theta }}$ denotes the vector form of the tuple $\left\{ {{{X,Y}},\gamma } \right\}$,and ${{{F}}_i}({\pmb{\theta }},\rho,{\pmb{\delta }} ) < 0$ represents the LMI constraints arising from Theorem 3. Note that while both $\rho $ and ${\pmb{\delta }} $ enter the RCP problem in the same manner,they act differently on the system; $\rho $ enters the control law as a scheduling parameter,while ${\pmb{\delta }} $ is treated as parameter perturbation whose adverse impact should be restrained by the control law. Therefore,(18) is indeed a robust optimal LPV control problem,taking both demands of gain scheduling and robustness into account.

To make the ISA algorithms applicable to problem (18),firstly it is necessary to define an equivalent scalar function as well as the associated sub-gradient. Let

where ${{{\tilde F}}_i}$ and ${{{\hat F}}_i}$$'$s dependence on ${\pmb{\theta }}$,$\rho $ and ${\pmb{\delta }} $ is omitted for brevity. Then a scalar feasibility indicator function^[22] is adopted. Define

where ${{{\hat F}}_ + }$ denotes the projection of matrix ${{\hat F}}$ onto the cone of positive semidefinite matrices,and the subscript ``fro'' denotes the Frobenius norm that we adopt in this paper. It is easy to verify that $\varphi ({\pmb{\theta }},\rho, {\pmb{\delta }} ) > {\rm{0}}$ if and only if ${{\hat F}}({\pmb{\theta }},\rho,{\pmb{\delta }} ) \le {\rm{0}}$ does not hold,and it is zero,otherwise. Define

Then,the sub-gradient of $\varphi ({\pmb{\theta }},{\pmb{\delta }} )$ is given as

Then we cope with the parameter dependency problem. Note that $\varphi ({\pmb{\theta }},\rho,{\pmb{\delta }} ) \le {\rm{0}}$ ( $\varphi ({\pmb{\theta }},\rho,{\pmb{\delta }} ) = {\rm{0}}$ in our case) must hold for all $\rho \in {\cal P}$. To handle the LMI$'$s dependence on $\rho $,a so-called gridding technique is commonly used[1,23]. After $L$ gridding points are selected,i.e., ${{\cal P}_g} = \left\{ {{\rho ^{(l)}}} \right\}_{l = 1}^{L} \subset {\cal P}$ and scenarios of ${\pmb{\delta }} $ are collected, problem (18) is further approximated as

which has totally $L \times N$ constraints.

In this paper,we use a different way to deal with the scheduling parameter,that we cast it as a role of the scenario,generating random samples rather than making grids upon it. This means the scheduling parameters act like uncertainty parameters during the computation,which can be easily understood by observing that $\rho $ and ${\pmb{\delta }} $ enter problem (18) in the same manner. So problem (18) is finally translated into the following scenario design problem that can be solved by ISA:

The benefit of doing so is obvious that randomly sampling scheduling parameter does not increase the sample size of scenarios,while greatly reducing the number of LMI constraints for each scenario, compared to the gridding approach. And if the solution of (19) exists,it is supposed to be better generalized to all $\rho \in {\cal P}$ than that of using the gridding approach,as it has experienced more instances of $\rho $.

In this section,the ISA proposed in the above section is applied to design an LPV longitudinal controller of a small helicopter for its cruise flight. We employ a Hirobo 90-class model scale helicopter,shown in Fig.1,as the testbed of our design method. A small helicopter is accurately modeled by the fuselage$'$s 6-DOF rigid body model augmented with the first-order main rotor and stabilizer bar$'$s flapping dynamics,whose complete model description includes 13 states: 3 linear velocity components in body-axis ($u,v,w$),3 angular rates ($p,q,r$),3 Euler angles ($\phi,\theta,\varphi $),2 flapping angles of main rotor$'$s tip-path-plane (TPP) (${a_1},{b_1}$) and 2 flapping angles of Bell-Hiller bar$'$s TPP ($c,d$). A minimum complexity nonlinear model^[24] is used both for linear model extraction and simulation verification. Earlier work on modeling a small helicopter with stabilizer bar^[25] showed that the main rotor and stabilizer bar$'$s flapping dynamics could be lumped and modeled just by main rotor$'$s 2-order TPP flapping dynamics,with the stabilizer bar$'$s effect counted into main rotor$'$s equivalent flapping time constant. So for longitudinal motion,the state vector is defined to be ${[\begin{array}{*{20}{c}} u&{{a_1}}&w&q&\theta \end{array}]^{\rm{T}}}$,where $u$ is the body axis forward speed,${a_1}$ is the longitudinal flapping angle,$w$ is the vertical speed,$q$ is the pitch rate and $\theta $ is the pitch angle.

Fig. 1

	Download: JPG larger image
Fig. 1. The equipped 90-class model helicopter in flight.

In order to describe the flight dynamics,many physical parameters need to be evaluated,some of which are related to aerodynamics,and rotor flapping are hard to be determined,as listed in Table I.

Table I

Table I UNCERTAIN PARAMETERS OF A SMALL HELICOPTER AND THEIR DISTRIBUTIONS Parameter Description Nominal value Distribution ${a_{{\rm{mr}}}}$ Main blade lift curve slope 5.5 rad$^{-1}$ U(5.0, 6.0) $C_{D0}^{{\rm{mr}}}$ Main blade zero lift drag coefficient 0.024 U(0.018, 0.03) ${\tau _{{\rm{mr}}}}$ Main rotor flapping time constant 0.1s U(0.06, 0.14) ${K_\beta }$ Hub torsion stiffness 55 N$\cdot$m/rad U(50, 70)

Table I UNCERTAIN PARAMETERS OF A SMALL HELICOPTER AND THEIR DISTRIBUTIONS

The listed items are the uncertain parameters we care about for controller design,whose nominal values are determined by citing the parameter data of an X-Cell 60-class helicopter^[25],and we assume they are uniformly distributed around the nominal values.

The model dynamically varies depending on the advanced ratio $\mu $,defined by $\mu = {u /{\Omega R}}$,where $\Omega $ and $R$ are main rotor$'$s revolution,speed and radius,respectively. It is reasonable to select $u$ to be the scheduling parameter, i.e.,$\rho = u$,as $\Omega $ is usually kept constant by an independent engine controller during the flight. The flight envelope of interest covers the speed regime of $0\sim25$m/s,so the scheduling parameter $\rho $ is assumed to be uniformly distributed in this interval.

Following the trim and linearization procedures,the helicopters$'$ longitudinal LPV model can be extracted from its nonlinear simulation model,expressed as (20) in the next page, where ${{\pmb{\delta }} _{{\rm{lon}}}}$ and ${{\pmb{\delta }} _{{\rm{col}}}}$ are longitudinal cyclic input,and collective input,respectively; ${\theta _{\rm{e}}}$ and ${\phi _{\rm{e}}}$ are trimmed angles; $X,D,Z,M$ are the stability derivatives and $K$ is the control derivative; and $y$ denotes the output vector.

Similar to the synthesis procedures for LTI systems,weighting functions are defined here to characterize the overall closed-loop performance. The open-loop interconnection of the corresponding weighted system is shown in Fig. 2,where

Fig. 2

	Download: JPG larger image
Fig. 2. Weighted open loop interconnection of helicopter plant.

From the interconnection shown in Fig. 2,it is easy to construct a 12-order augmented open-loop LPV model of the form (13) for the longitudinal dynamics. The corresponding SCP has 157 design variables. To get an initial guess of the design variables,standard SA is implemented on 1000 scenarios using YALMIP^[26] under Matlab environment,which nearly uses up the computational capacity of a personal computer equipped with 2GB memory. The calculated optimal performance level is 2.09,and the initial guess of $\gamma $ is set to be 6.0,allowing some gap above the calculated value. The initial hypercube constraints,centered on the initial guess mentioned above,are chosen to be large enough to guarantee the confidence that the optima are included in the hypercube,whose volume is about $2.62 \times {10^{365}}$.

We set $\tau = 0.07$,$\varepsilon = 0.05$,and $\beta = {10^{ - 5}}$ to execute the proposed two ISA algorithms. According to (3),6742 scenarios are created. For Algorithm 2, we select $\eta = 0.95$,

and $K = 400$,and then using (8),the repeat times of sampling inspection,$M$,is set to be 51. For comparison,we also implement in the same settings,the SE algorithm based on the MVE cutting-plane method,which is proposed by Wada et al. ^[18] and regarded as the most sophisticated framework among current SE methods. All the algorithms are implemented on a personal computer with 2.26GHz CPU and 2GB memory. Table II shows the process information of the three algorithms.

Table Ⅱ

Table Ⅱ PROCESS INFORMATION OF DIFFERENT ALGORITHMS Algorithm Scenarios used per iteration Number of iterations Time(s) Total amount of scenarios Optimum(${\gamma ^ * }$ ) 1 6742 595 28955 (7416) 6742 3.5229 2 400 512 11036 (7416) 6742 3.1323 Wada$'$s 397 423 24517 (21705) 19732 3.3242

Table Ⅱ PROCESS INFORMATION OF DIFFERENT ALGORITHMS

In Table II,the execution time (the 4th column) of the algorithms involves both the time spent on the iterations and the time to generate enough scenarios for the algorithm; the latter is particularly listed in the brackets to facilitate the comparison. As observed from the table,Algorithm 1 is the least efficient in the respect of execution time; the scenarios consumed by Wada$'$s algorithm are nearly 3 times that of Algorithms 1 and 2,thus making it rather sluggish; the batch algorithm (Algorithm 2) spent just 45$\%$ of the time taken by Wada$'$s algorithm.

We apply Monte Carlo approach to verify the robust feasibility of the solutions obtained from different algorithms. Firstly,we estimate the empirical probability of robust feasibility for both nominal solution and the solutions by the three algorithms mentioned above,by checking whether the LMIs (16a)$\sim$(16c) (or $\varphi ({\pmb \delta},\rho,{{\pmb \theta }}) \le 0$ ) hold for each test scenario. Secondly,to verify the actual closed-loop robustness,we adopt the closed-loop LMI formula^[20] that conditions the $\gamma $-level L2-norm performance,given by

where

and ${{W}}$ is the matrix associated with the Lyapunov function of the closed-loop system,defined here by

Ten thousand scenarios are generated for the verifications,and the results are shown in Table III. For comparison purpose,we also put into the table the results obtained from other three algorithms: LPV control (LPV-C) on nominal plant,SA on 1000 scenarios (whose solution is used as the initial guess of the iterative algorithms), and Wada$'$s algorithm.

Table Ⅲ

Table Ⅲ FEASIBILITY VERIFICATION RESULTS OF DIFFERENT ALGORITHMS Solution number Algorithm $\Pr \left( {{\pmb \delta} \in \Delta, \varphi ({\pmb \delta}, \rho, {{\pmb \theta }}) \le 0} \right)$$\left( \% \right)$ Closed-loop feasibility$\left( \% \right)$ 1 LPV-C on nominal 0.01 0 2 SA on 1\,000 scenarios 70 0 3 1 99.98 99.46 4 2 100 99.86 5 3 99.91 99.84

Table Ⅲ FEASIBILITY VERIFICATION RESULTS OF DIFFERENT ALGORITHMS

From Table III,it seems that the robustness improvement is encouraging after performing the first kind of verification,as the empirical estimates of probabilistic robust feasibility for both Algorithms 1 and 2 are much better than their theoretical bounds. However,the results of closed-loop LMI verification are a bit disappointing. Recall that in order to make the controller implementable,we construct it based on the nominal plant,see (17) where we use ${{A}}(\rho )$ rather than ${{A}}(\rho,{\pmb \delta})$. As discussed in Section III,both $\rho $ and ${\pmb \delta}$ act in the same way from a pure mathematical point of view,and in the general SQLF settings, ${{A}}( \cdot )$ should depend on both $\rho $ and ${\pmb \delta}$. So neglecting its dependence on the uncertainty ${\pmb \delta}$ results in a loss of robustness. But from a practical point of view,things are still acceptable. Similar results were reported in ^[23],but the performance degradation of ours is remarkably smaller.

Finally,the closed-loop performance of robust LPV controller constructed using the solution of Algorithm 2 is verified through simulations on helicopter$'$s nonlinear model. To construct the controller as formulated in (17),nominal LPV matrices ${{A}}(\rho )$ and ${{B}}(\rho )$ are approximated by polynomial matrices driven by $\rho $. We exert step-like reference signal to either forward or vertical channels respectively,while keeping the other at zero state. Ten uncertain nonlinear models are used for the verification,and the results are shown in Figs.3 and 4,where the dashed lines represent the desired responses (output of the reference model) and the solid lines represent the actual responses. Note that the nominal LPV controller cannot pass the closed-loop feasibility verification for any scenario candidate,as shown in Table III,so there is no need to verify it on the nonlinear model again.

From Figs.3 and 4 seen that in the presence of model perturbations,the closed-loop system tracks reference model outputs very well. This means that governed by the LPV controller, the system consistently presents a second-order closed-loop characteristic specified by the reference model in a wide speed range (0$\sim$23m/s),which is quite desirable for the high-layer motion planning design. It is further observed from Fig.3 that the variation of vertical speed is kept very small during the high-speed forward flight,which means that the longitudinal phugoid mode is eliminated,demonstrating that the decoupling performance of the controller is satisfactory.

Fig. 3

	Download: JPG larger image
Fig. 3. Forward and vertical speed responses to step-like excitation on the forward speed channel.

Fig. 4

	Download: JPG larger image
Fig. 4. Forward and vertical speed responses to step-like excitation on the vertical speed channel.

For a complex system like a small helicopter,neither trim nor linearization procedures are simple,as seen from Table I it takes more than 1s to create a scenario. In this situation,those SE algorithms are too time-consuming. The proposed batch algorithm (Algorithm 2) adopts both advantages of Algorithm 1 and SE algorithms, i.e.,it consumes the same amount of scenarios in all as the former,and queries as many as the latter for feasibility verifications.

Numerical stability is a problem of optimization algorithms using cutting-plane techniques. For the controller design of a small helicopter$'$s longitudinal motion,a 157-dimensional optimization problem,numerical stability would gradually deteriorate,i.e., the condition numbers of the matrices generated from the algorithm become much larger,as the cutting-plane-based algorithm proceeds. The algorithm might stop unexpectedly with fatal errors,if we cancel the constraint pruning codes. The aforementioned heuristic logic for removing redundant constraints can effectively improve the numerical stability as no fatal error occurs once the logic works. However,the curse of dimension is still a problem that we cannot thoroughly circumvent. For SCP problems of larger-scale, e.g.,larger than 157-dimensions,a more elaborate constraint pruning technique would be necessary.

In this paper,we have presented iterative algorithms for the SA. For the robust control of a complex LPV system which is time-costly to create through linearization,the proposed batch algorithm has been shown to be more efficient than the popular SE methods. In the robust LPV control context,we have also proposed to cope with parameter dependency using the randomization technique identically on both uncertainty and scheduling parameters,and to prune the redundant constraints using a heuristic logic based on relevance ratio ranking. All these techniques together comprise a complete method,which has been applied to design robust LPV controller for small helicopter$'$s cruise flight. Multiple uncertain physical parameters have been considered. Simulation studies have been provided to demonstrate its effectiveness.