2. Department of Control Science and Engineering, Zhejiang University, Hangzhou 310027, China;
3. School of Aeronautics and Astronautics, Zhejiang University, Hangzhou 310027, China and Department of Control Science and Engineering, Zhejiang University, Hangzhou 310027, China
Gain scheduling based on linear parametervarying (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 socalled scheduling parameters that can be measured online or apriori known. Then based on these LPV models,the controller is synthesized in a linearlike 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 Lyapunovbased 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 timeconsuming.
SA,developed by Calafiore and Campi^{[11, 12]},gives a oneshot 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 firstorder polynomial to the number of design variables^{[13]}. So when facing highdimensional 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 threefold. 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 commonlyused 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 outputfeedback 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 outputfeedback 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.
Ⅱ. ISA A. Standard SAA general robust convex program (RCP) takes the following form:
$ \begin{array}{l} {\rm{RCP}}:\mathop {\;\min }\limits_{{{\pmb \theta }} \in {{\Theta }}} \;{{\pmb c}^{\rm{T}}}{\pmb{\theta }},\;\;\;\;\mbox{s. t.}\quad \varphi ({\pmb{\theta }},{\pmb{\delta }} ) \le 0,\;\forall {\pmb{\delta }} \in {{\Delta }},\label{eq:rel1} \end{array} $  (1) 
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 NPhard problem due to the possible infinity of constraints.
SA looses the constraint satisfaction from ``all" to ``almost" by applying the socalled 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
$ \begin{array}{l} {\rm{SCP}}_N:\;\;\mathop {\min }\limits_{{\pmb{\theta }} \in {{\Theta }} } \;{{\pmb c}^{\rm{T}}} {\pmb{\theta }},\;\;\;\mbox{s. t.}\quad \varphi ({\pmb{\theta }},{{\pmb{\delta }} ^{(k)}}) \le 0,\;k = 1,\cdots,N. \end{array} $  (2) 
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
$ N \ge \frac{2}{\varepsilon }\left( {\log \frac{1}{\beta } + {n_{\pmb{\theta }}}} \right), $  (3) 
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 oneshot 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.
B. ISAInspired 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
$E(Q,{{\theta }_{c}})=\{\theta \in {{\mathbf{R}}^{{{n}_{\theta }}}}:\theta =Qx+{{\theta }_{c}},\left\ x \right\<1,x\in {{\mathbf{R}}^{{{n}_{\theta }}}}\},$ 
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
${E_{{\rm{mv}}}}({{P}}) = \mathop {\arg }\limits_{Q,{{\pmb{\theta }}_{{c}}}} \max \{ {\rm{vol}}(E({ Q},{{\pmb{\theta }}_{{c}}})),E({Q},{{\pmb{\theta }}_{{c}}}) \subset {{P}}\}$. 
It is easily seen that the MVE can be found by solving the following optimization problem:
$\begin{array}{l} \max \;\ln \left( {\det {Q}} \right),~~~\mbox{s. t.}\,\left\ {{\pmb a}_i^{\rm{T}}{Q}} \right\ + {\pmb a}_i^{\rm{T}}{{\pmb{\theta }}_{{c}}} \le {{b}_i},\;\;i = {\rm 1}, \cdots,L. \end{array}$ 
To facilitate the planecutting 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
$ {{{\Theta }} _{\rm 1}}: = \{ {\pmb{\theta }} \in {{\bf{R}}^{{n_{\pmb{\theta }}}}}:  {{ r}_i} \le {\pmb{e}}_i^{\rm{T}}({\pmb{\theta }}  {{\pmb{\theta }}_{\rm 1}}) \le {{r}_i},i = {\rm 1},\cdots,{n_{\pmb{\theta }}}\}, $  (4) 
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
$ {{P}_1}: = \left\{ {{\pmb{\theta }} \in {{\bf{R}}^{{n_{\pmb{\theta }}}}}:{{{A}}_1}{\pmb{\theta }} \le {{{B}}_1}} \right\}. $ 
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
$ j = j + 1,\;\;{{{\hat {\pmb{\theta }} }}_j} = {{\pmb{\theta }}_l},\; $ 
$ f = 0,\;\;{\pmb {g}} = {\pmb {c}}, $ 
else set
$ f = \varphi ({{\pmb{\theta }}_l},{{\pmb{\delta }} ^{(k)}}),\;{\pmb {g}} = {\partial _{\pmb{\theta }}}\varphi ({{\pmb{\theta }}_l},{\pmb{\theta }} ). $ 
2) Update. Add a new cutting plane to ${P_l}$ :
$ {{{A}}_{l + 1}} = \left[{\begin{array}{*{20}{c}} {{{{A}}_l}}\\ {{{\pmb{g}}^{\rm{T}}}} \end{array}} \right],\;\;{{\pmb{B}}_{l + 1}} = \left[{\begin{array}{*{20}{c}} {{{{B}}_l}}\\ {{{\pmb{g}}^{\rm{T}}}{{\pmb{\theta }}_l}  f} \end{array}} \right],\;\;\;\; $  (5) 
find the MVE ${E_{l + 1}}$ in the new polytope:
$ \;{{P}_{l + 1}}: = \left\{ {{\pmb{\theta }} \in {{\bf{R}}^{{n_{\pmb{\theta }}}}}: {{{A}}_{l + 1}}{\pmb{\theta }} \le {{{B}}_{l + 1}}} \right\}, $ 
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 constraintpruning 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
$ { R}({{\pmb{a}}_i},{{ b}_i}) = \frac{{{b_i}  {\pmb{a}}_i^{\rm{T}}{{\pmb{\theta }}_{{c}}}}}{{{{\left\ {{{{F}}^{\rm{T}}}{{\pmb{a}}_i}} \right\}_2}}}, $  (6) 
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
$ \begin{array}{l} {\rm{if}}\;l \le {n_{\pmb{\theta }}},\\ \qquad {{{A}}_{l + 1}} = \left[{\begin{array}{*{20}{c}} {{{{A}}_l}}\\ {{{\pmb{g}}^{\rm{T}}}} \end{array}} \right],\;{{{B}}_{l + 1}} = \left[{\begin{array}{*{20}{c}} {{{{B}}_l}}\\ {{{\pmb{g}}^{\rm{T}}}{{\pmb{\theta }}_l}  f} \end{array}} \right];\\ {\rm{else}},\;\\ \qquad {{{A}}_{l + 1}} = {{{A}}_l},\;{{{B}}_{l + 1}} = {{{B}}_l},\\ \;\;\;\;\;\;\;\;\;n = \mathop {\arg }\limits_i \max \{R({{\pmb{a}}_i},{b_i}),\;i = 1,2,\cdots,3{n_{\pmb{\theta }}}\} ,\;\\ \qquad {{\pmb{a}}_n} = {{\pmb{g}}^{\rm{T}}},\;{b_n} = {{\pmb{g}}^{\rm{T}}}{{\pmb{\theta }}_l}  f,{\rm{ }} \end{array} $  (7) 
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 MVEbased 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 suboptimal 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 randomlysampled 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
$ \eta \left( { {{N},k,{m}}} \right) = \frac{1}{{\sum\limits_{i = k}^{N} {{{ ({\frac {i}{{N}}} )}^ {{m}k}}} }}. $  (8) 
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.,
$ l\left( {x\lambda } \right) = \left( {\begin{array}{*{20}{c}} k\\ x \end{array}} \right){\lambda ^x}{(1  \lambda )^{k  x}},\; x = 0,1,\cdots,k, $  (9) 
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
$ l(\{ {k^{(1)}},{k^{(2)}},\cdots,{k^{(m)}}\} \lambda ) = {\lambda ^{mk}}. $  (10) 
In the Bayesian sense,the posterior probability of $\lambda $ is given by
$ \Pr (\lambda \{ {k^{(\rm 1)}},{k^{(\rm 2)}},\cdots,{k^{(m)}}\} ) =\notag\\ \quad\frac{{l(\{ {k^{(\rm 1)}},{k^{(\rm 2)}},\cdots,{k^{(m)}}\} \lambda )g(\lambda )}}{{\sum\limits_\lambda {l(\{ {k^{(\rm 1)}},{k^{(\rm 2)}},\cdots,{k^{(m)}}\} \lambda )g(\lambda )} }}, $  (11) 
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
$ g(\lambda ) = \left\{ {\begin{array}{*{20}{l}} \frac{1}{{{N} + 1  k}},\\ {0}, \end{array}\begin{array}{*{20}{l}} {\lambda \in \left\{ {\frac{k}{N},\frac{{k + 1}}{N},\cdots,1} \right\},}\\ {{\mbox{otherwise}}}. \end{array}} \right. $  (12) 
Inserting (10) and (12) into (11),the confidence that $\lambda $ equals to 1 is then obtained as
$ \eta \left( {{N},k,m} \right) = \Pr (\lambda = 1\{ {k^{(1)}},{k^{(\rm 2)}},\cdots,{k^{(m)}}\} )=\\ \quad\quad\frac{1}{{\sum\limits_{i = k}^{N} {{{( {{{\frac {i}{N}}}} )}^{{m}k}}}. }} $ 
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
$ j = j + 1,\;\;{{\bf{\hat {\pmb{\theta }} }}_j} = {{\pmb{\theta }}_l},\; $ 
$ f = 0,\;\;{\pmb{g}} = {\pmb{c}}; $ 
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:
$ \;{{ P}_{l + 1}}: = \left\{ {{\pmb{\theta }} \in {{\bf{R}}^{{n_{\pmb{\theta }}}}}:{{{A}}_{l + 1}}{\pmb{\theta }} \le {{{B}}_{l + 1}}} \right\}, $ 
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$ randomlypicked $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.
Ⅲ. ROBUST OUTPUTFEEDBACK CONTROL OF LPV SYSTEM USING ISA A. Preliminary on Robust LPV ControlGiven 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.,
$ \;\left[{\begin{array}{*{20}{c}} {{{\dot x}}}\\ {{z}}\\ {{y}} \end{array}} \right] = \left[{\begin{array}{*{20}{c}} {{{A}}\left( {\rho,{\pmb{\delta }} } \right)} &{{{{B}}_{\bf{1}}}\left( {\rho,{\pmb{\delta }} } \right)}&{{{{B}}_{\bf{2}}}\left( {\rho,{\pmb{\delta }} } \right)}\\ {{{{C}}_{\bf{1}}}\left( {\rho,{\pmb{\delta }} } \right)}&{{{{D}}_{{\bf{11}}}}\left( {\rho,{\pmb{\delta }} } \right)} &{{{{D}}_{{\bf{12}}}}\left( {\rho,{\pmb{\delta }} } \right)}\\ {{{{C}}_{\bf{2}}}\left( {\rho,{\pmb{\delta }} } \right)}&{{{{D}}_{{\bf{21}}}}\left( {\rho,{\pmb{\delta }} } \right)}&{{{{D}}_{{\bf{22}}}}\left( {\rho,{\pmb{\delta }} } \right)} \end{array}} \right]\;\left[{\begin{array}{*{20}{c}} {\pmb{x}}\\ {\pmb{d}}\\ {\pmb{u}} \end{array}} \right], $  (13) 
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 apriori 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
$ \left[{\begin{array}{*{20}{c}} {{{{\pmb{\dot x}}}_{{k}}}}\\ {\pmb{u}} \end{array}} \right] = \left[{\begin{array}{*{20}{c}} {{{{A}}_{{k}}}\left( \rho \right)}&{{{{B}}_{{k}}}\left( \rho \right)}\\ {{{{C}}_{{k}}}\left( \rho \right)}&{{{{D}}_{{k}}}\left( \rho \right)} \end{array}} \right]\left[{\begin{array}{*{20}{c}} {{{\pmb{x}}_{{k}}}}\\ {\pmb{y}} \end{array}} \right]. $  (14) 
In this paper,we are concerned with the induced L2norm performance of closedloop LPV system (13) and (14),which is defined as
$ {\left\ {{{{G}}_{{F_{\cal P}}}}} \right\_{{\rm{i}},2}} = \;\mathop {\sup }\limits_{\rho \in {F_{\cal P}}} \mathop {\sup }\limits_{{{\left\ {{\pmb d}} \right\}_2} \ne 0,{{\pmb d}} \in {{{L}}^{\bf{2}}}} \frac{{{{\left\ {{\pmb z}} \right\}_2}}}{{{{\left\ {{\pmb d}} \right\}_2}}}, $ 
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 loopshifting 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 normpreserving transformations on ${\pmb z}/{\pmb d}$ and invertible transformations on ${\pmb y}/{\pmb u}$:
$ \begin{align} &\left[{\begin{array}{*{20}{c}} {{\pmb{\dot x}}}\\ {{{\pmb{z}}_{\bf{1}}}}\\ {{{\pmb{z}}_{\bf{2}}}}\\ {\pmb{y}} \end{array}} \right] = \notag\\ &~\left[{\begin{array}{*{20}{c}} {{{A}}\left( {\rho,{\pmb{\delta }} } \right)}&{{{{B}}_{11}}\left( {\rho,{\pmb{\delta }} } \right)}&{{{{B}}_{12}}\left( {\rho,{\pmb{\delta }} } \right)} &{{{{B}}_2}\left( {\rho,{\pmb{\delta }} } \right)}\\ {{{{C}}_{11}}\left( {\rho,{\pmb{\delta }} } \right)}&{\bf{0}}&{\bf{0}}&{\bf{0}}\\ {{{{C}}_{12}}\left( {\rho,{\pmb{\delta }} } \right)}&{\bf{0}}&{\bf{0}}&{{{\pmb {I}}_{{n_{\pmb u}}}}}\\ {{{{C}}_2}\left( {\rho,{\pmb{\delta }} } \right)}&{{\bf 0}}&{{{\pmb{I}}_{{n_{\pmb y}}}}}&{\bf{0}} \end{array}} \right]\left[{\begin{array}{*{20}{c}} {\pmb{x}}\\ {{{\pmb{d}}_{\bf{1}}}}\\ {{{\pmb{d}}_{\bf{2}}}}\\ {\pmb{u}} \end{array}} \right]. \end{align} $  (15) 
For the L2norm robust control of an LPV system in the presence of model uncertainty,we give a sufficient condition as follows
. Theorem 3. Given the uncertain LPV system$\sum\nolimits_{{\cal P},{{\Delta }} }^S {} $,and a controller ${K_{\cal P}}$,the closedloop system is guaranteed to have an L2norm 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:
$ \left[\begin{array}{*{20}{cccc}} E_{\pmb x}(\rho,{\pmb \delta } )&{XC}_{11}^{\rm T}(\rho,{\pmb \delta} )&{B _1}(\rho, {\pmb \delta } )\\ C_{11}(\rho,{\pmb \delta} )X&\gamma {\pmb I}_{n_{z_1}}&{\bf 0}\\ {B_1^{\rm T}(\rho,{\pmb \delta } )}&{\bf 0}&{  \gamma {\pmb{I}}_{{n}_{\pmb d}}} \end{array} \right] < 0, $  (16a) 
$ \left[{\begin{array}{*{20}{cccc}} {{{{E}}_{\bf{y}}}(\rho, {\pmb{\delta }} )}&{{{Y}}{{{B}}_{11}}(\rho,{\pmb{\delta }} )} &{{{C}}_{\rm{1}}^{\rm{T}}(\rho,{\pmb{\delta }} )}\\ {{{B}}_{{\rm{11}}}^{\rm{T}}(\rho,{\pmb{\delta }} ){{Y}}}&{  \gamma {{\pmb{I}}_{{n_{{{\pmb d}_1}}}}}}&{\bf{0}}\\ {{{{C}}_1}(\rho,{\pmb{\delta }} )}&{\bf{0}}&{  \gamma {{\pmb{I}}_{{n_{\pmb z}}}}} \end{array}} \right] < 0, $  (16b) 
$ \left[{\begin{array}{*{20}{ccc}} {{X}}&{{{\pmb{I}}_n}}\\ {{{\pmb{I}}_n}}&{{Y}} \end{array}} \right] \ge 0, $  (16c) 
where
$ {{{E}}_{{\pmb x}}}(\rho,{\pmb{\delta }} ) = {{X}}{{{\hat A}}^{\rm{T}}} (\rho,{\pmb{\delta }} ) + {{\hat A}}(\rho,{\pmb{\delta }} ){{X}}  \gamma {{{B}}_2}(\rho,{\pmb{\delta }} ){{B}}_2^{\rm{T}}(\rho, {\pmb{\delta }} ), $ 
$ {{{E}}_{{\pmb y}}}(\rho,{\pmb{\delta }} ) = {{{\tilde A}}^{\rm{T}}}(\rho,{\pmb{\delta }} ) {{Y}} + {{Y\tilde A}}(\rho,{\pmb{\delta }} )  \gamma {{C}}_2^{\rm{T}}(\rho,{\pmb{\delta }} ){{{C}}_2}(\rho, {\pmb{\delta }} ), $ 
$ {{\hat A}}(\rho,{\pmb{\delta }} ) = {{A}}(\rho,{\pmb{\delta }} )  {{{B}}_2}(\rho,{\pmb{\delta }} ){{{C}}_{12}}(\rho,{\pmb{\delta }} ), $ 
and
$ {{\tilde A}}(\rho,{\pmb{\delta }} ) = {{A}}(\rho,{\pmb{\delta }} )  {{{B}}_{12}}(\rho,{\pmb{\delta }} ) {{{C}}_2}(\rho,{\pmb{\delta }} ). $ 
The theorem can be easily proved by directly applying Wu$'$s theorem (see ^{[20]},Theorem 3.3.1) on parameter dependent outputfeedback 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 parameterdependent Lyapunov function (PDLF). Once such ${{X}}$ and ${{Y}}$ are found,we can then construct the robust outputfeedback LPV controller based on the nominal statespace matrices. Define^{[20]}
$ \begin{array}{c} {{Q}}(\rho ) = {\gamma ^{  1}}\left( {{{Y}}  {{{X}}^{  1}}} \right) > 0,\\ \\ {{F}}(\rho ) =  [\gamma {{B}}_2^{\rm T}(\rho ){{{X}}^{  1}} + {{D}}_{12}^{\rm T}{{{C}}_1}(\rho )],\\ \\ {{L}}(\rho ) =  [\gamma {{{Y}}^{  1}}{{C}}_2^{\rm T}(\rho ) + {{{B}}_1}(\rho ){{D}}_{21}^{\rm T}],\\ \\ {{H}}(\rho ) =  [\gamma {{{X}}^{  1}}{{{A}}_F}(\rho ) + \gamma {{A}}_F^{\rm T}(\rho ){{{X}}^{  1}} + {{C}}_F^{\rm T}(\rho )\times\\ \\ {{{C}}_F}(\rho ) + {{{X}}^{  1}}{{{B}}_1}(\rho ){{B}}_1^{\rm T}(\rho ){{{X}}^{  1}}], \end{array} $ 
where
$ {{{A}}_F}(\rho ) = {{A}}(\rho ) + {{{B}}_2}(\rho ){{F}}(\rho ) $ 
and
$ {{{C}}_F}(\rho ) = {{{C}}_1}(\rho ) + {{{D}}_{12}}{{F}}(\rho ). $ 
And further,let
$ \begin{align*}&{{M}}(\rho ) = {{H}}(\rho ) + \\ &\quad{\gamma ^2}{{Q}}(\rho )[ {{{Q}}^{  1}}(\rho ){{YL}}(\rho ){{{D}}_{21}}  \gamma {{{B}}_1}(\rho )]{{B}}_1^{\rm T}(\rho ){{{X}}^{  1}},\end{align*} $ 
then the $n$order strictly proper controller ${K_{\cal P}}$ achieving $\gamma$performance is presented as
$ \begin{align} & {{A}_{k}}(\rho )=A(\rho )+{{B}_{2}}(\rho )F(\rho )+{{\gamma }^{1}}{{Q}^{1}}(\rho )YL(\rho ){{C}_{2}}(\rho )~~ \\ & ~{{\gamma }^{2}}{{Q}^{1}}(\rho )M(\rho ),\\ & {{B}_{k}}(\rho )={{\gamma }^{1}}{{Q}^{1}}(\rho )YL(\rho ),\\ \end{align}$  (17) 
${{{C}}_{{k}}}(\rho ) = {{F}}(\rho ),$ 
${{{D}}_{{k}}}(\rho ) = {\pmb{0}}.$ 
The RCP problem associated with Theorem 3 can be expressed as
$ \begin{array}{l} \mathop {\min }\limits_{{\pmb{\theta }} \in {{\Theta }} } \;\gamma \;\;\;\;\mbox{s. t.} \;{{{ F}}_i}({\pmb{\theta }},\rho, {\pmb{\delta }} ) < 0,\;i = {\rm{1}},{\rm{2}},{\rm{3}},\;\forall \rho \in { {\cal P}},{\pmb{\delta }} \in {{\Delta }},\end{array} $  (18) 
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.
C. Computational Issues for Applying ISATo make the ISA algorithms applicable to problem (18),firstly it is necessary to define an equivalent scalar function as well as the associated subgradient. Let
$ {{\hat F}}({\pmb{\theta }},\rho,{\pmb{\delta }} ) = {\rm{diag}}\{{{{\tilde F}}_{\rm{1}}},{{{\tilde F}}_{\rm{2}}}, \cdots,{{{\tilde F}}_m}\} = {{{\hat F}}_{\rm{0}}} + \sum\limits_{i = {\rm{1}}}^{{n_{\pmb{\theta }}}} {{{\pmb{\theta }}_i}{{{{\hat F}}}_i}}, $ 
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
$ \varphi ({\pmb{\theta }},\rho,{\pmb{\delta }} ): = {\left\ {\hat F}_+(\pmb{\theta },\rho ,\pmb{\delta }) \right\_{\rm fro}}, $ 
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
$ \nabla \varphi ({\pmb{\theta }},{\pmb{\delta }} ): = \frac{{\rm{1}}}{{\varphi ({\pmb{\theta }}, {\pmb{\delta }} )}}\left[{\begin{array}{*{20}{c}} {{\rm{tr(}}{{{F}}_1}{{{F}}_ + }({\pmb{\theta }},{\pmb{\delta }} ))}\\ \cdots \\ {{\rm{tr(}}{{{F}}_{{n_{\pmb{\theta }}}}}{{{F}}_ + }({\pmb{\theta }},{\pmb{\delta }} ))} \end{array}} \right]. $ 
Then,the subgradient of $\varphi ({\pmb{\theta }},{\pmb{\delta }} )$ is given as
$ {\partial _{\pmb{\theta }}}\varphi ({\pmb{\theta }},{\pmb{\delta }} ) = \left\{ {\begin{array}{*{20}{c}} {\nabla \varphi ({\pmb{\theta }},{\pmb{\delta }} ),\;\;{\rm{if}}\;\varphi ({\pmb{\theta }},{\pmb{\delta }} ) \ne 0},\\ {0,\;\;\;\;\;\;\;\;\;\;\; {\rm{otherwise.}}\;\;\;\;\;} \end{array}} \right. $ 
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 socalled 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
$ \begin{array}{l} \mathop {\min }\limits_{{\pmb{\theta }} \in \Theta } \;\gamma,\;\;\;\;\mbox{s. t.}\quad\varphi ({\pmb{\theta }},{\rho ^{(l)}},{{\pmb{\delta }} ^{(n)}}) \le 0,\;\\ \;\;\;\;\;\;l = 1,\cdots,{L},\;n =1,\cdots,N, \end{array} $ 
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:
$ %\begin{array}{llll} {\rm{SCP}}_N:\mathop {\rm min}\limits_{{\pmb{\theta }} \in {{\Theta }} } \;\gamma,\;\; \mbox{s. t.}\;\; %\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \varphi ({\pmb{\theta }},{\rho^{(n)}},{{\pmb{\delta }} ^{(n)}}) \le 0,\;\;\;n = 1,\cdots,N. %\end{array} $  (19) 
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 $.
Ⅳ. APPLICATION TO THE CONTROL OF A SMALL HELICOPTER A. Model SetupIn 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 90class model scale helicopter,shown in Fig.1,as the testbed of our design method. A small helicopter is accurately modeled by the fuselage$'$s 6DOF rigid body model augmented with the firstorder main rotor and stabilizer bar$'$s flapping dynamics,whose complete model description includes 13 states: 3 linear velocity components in bodyaxis ($u,v,w$),3 angular rates ($p,q,r$),3 Euler angles ($\phi,\theta,\varphi $),2 flapping angles of main rotor$'$s tippathplane (TPP) (${a_1},{b_1}$) and 2 flapping angles of BellHiller 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 2order 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.
Download:


Fig. 1. The equipped 90class 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.
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 XCell 60class 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.
B. Controller DesignSimilar to the synthesis procedures for LTI systems,weighting functions are defined here to characterize the overall closedloop performance. The openloop interconnection of the corresponding weighted system is shown in Fig. 2,where
Download:


Fig. 2. Weighted open loop interconnection of helicopter plant. 
$ \begin{array}{l} {{{M}}_{{\rm{ref}}}} = {\rm{diag}}\left\{ {\begin{array}{*{20}{c}} {\frac{{(  0.05s + 1)}}{{{{{s}}^{\rm{2}}} + 1.6s + 1}}},&{\frac{1}{{s + 1}}} \end{array}} \right\},\\ \\ {{{W}}_{\pmb{\delta }} } = {\rm{diag}}\left\{ {\begin{array}{*{20}{c}} {\frac{{0.5s + 17.3}}{{s + 0.0577}}},&{\frac{{0.5s + 17.3}}{{s + 0.0577}}} \end{array}} \right\},\\ \\ {{{W}}_{{\pmb z}}} = {\rm{diag}}\left\{ {\begin{array}{*{20}{c}} {\frac{{10s + 100}}{{s + 1000}}},&{\frac{{10s + 100}}{{s + 1000}}} \end{array}} \right\},\\ \\ {{{W}}_{n1}} = {{{W}}_{n2}} = {\rm{diag}}\left\{ {\begin{array}{*{20}{c}} {0.0001},&{0.0001} \end{array}} \right\}. \end{array} $ 
From the interconnection shown in Fig. 2,it is easy to construct a 12order augmented openloop 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$,
$ \begin{align} & \left[\begin{matrix} {\dot{u}} \\ {{{\dot{a}}}_{1}} \\ {\dot{w}} \\ {\dot{q}} \\ {\dot{\theta }} \\ \end{matrix} \right]=\left[\begin{matrix} {{X}_{u}}(\rho ,\delta ) & g/g\cos {{\phi }_{\text{e}}} & 0 & {{u}_{\text{e}}}\tan {{\theta }_{\text{e}}} & g\cos {{\theta }_{\text{e}}} \\ {{D}_{u}}(\rho ,\delta ) & 1/1{{\tau }_{\text{mr}}}(\delta ) & {{D}_{w}}(\rho ,\delta ) & 1 & 0 \\ {{Z}_{u}}(\rho ,\delta ) & {{Z}_{{{a}_{1}}}}(\rho ,\delta ) & {{Z}_{w}}(\rho ,\delta ) & {{Z}_{q}}(\rho ,\delta ) & g\cos {{\phi }_{\text{e}}}\sin {{\theta }_{\text{e}}} \\ {{M}_{u}}(\rho ,\delta ) & {{M}_{{{a}_{1}}}}(\rho ,\delta ) & {{M}_{w}}(\rho ,\delta ) & {{M}_{q}}(\rho ,\delta ) & 0 \\ 0 & 0 & 0 & \cos {{\phi }_{\text{e}}} & 0 \\ \end{matrix} \right]\left[\begin{matrix} u \\ {{a}_{1}} \\ w \\ q \\ \theta \\ \end{matrix} \right]+\left[\begin{matrix} 0 & 0 \\ {{K}_{{{\delta }_{\text{lon}}}}}(\rho ,\delta ) & 0 \\ 0 & {{K}_{{{\delta }_{\text{col}}}}}(\rho ,\delta ) \\ 0 & K_{q}^{{{\delta }_{\text{col}}}}(\rho ,\delta ) \\ 0 & 0 \\ \end{matrix} \right] \\ & \times \left[\begin{matrix} {{\delta }_{\text{lon}}} \\ {{\delta }_{\text{col}}} \\ \end{matrix} \right] \\ & y={{\left[\begin{matrix} \ \ \ u & \ \ \ \ w\ \ \ \ & q\ \ \ \ & \theta \\ \end{matrix}\ \ \ \right]}^{\text{T}}}.n \\ \end{align} $  (20) 
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 cuttingplane 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.
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.
C. Simulation ResultsWe 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 closedloop robustness,we adopt the closedloop LMI formula^{[20]} that conditions the $\gamma $level L2norm performance,given by
$ \left[\!\!\! {\begin{array}{*{20}{c}} {{{A}}_{{\rm{clp}}}^{\rm{T}}(\rho,{\pmb \delta}){{W}} + {{W}}{{{A}}_{{\rm{clp}}}}(\rho,{\pmb \delta} )}&\!\!{{{W}}{{{B}}_{{\rm{clp}}}} (\rho,{\pmb \delta})}&{{{C}}_{{\rm{clp}}}^{\rm{T}}(\rho,{\pmb \delta})}\\ {{{B}}_{{\rm{clp}}}^{\rm{T}}(\rho,{\pmb \delta}){{W}}}&{  \gamma {{\pmb{I}}_{{n_{\pmb d}}}}} &{{{D}}_{{\rm{clp}}}^{\rm{T}}(\rho,{\pmb \delta})}\\ {{{{C}}_{{\rm{clp}}}}(\rho,{\pmb \delta})}&{{{{D}}_{{\rm{clp}}}}(\rho,{\pmb \delta})}&{  \gamma {{\pmb{I}}_{{n_{\pmb e}}}}} \end{array}}\!\!\! \right]\!<\!0,\notag $  (21) 
where
$ {{{A}}_{{\rm{clp}}}}(\rho,{\pmb \delta}) = \left[ {\begin{array}{*{20}{c}} {{{A}}(\rho,{\pmb \delta})}&{{{{B}}_2}(\rho,{\pmb \delta}){{{C}}_{{k}}}(\rho )}\\ {{{{B}}_{{k}}}(\rho,{\pmb \delta}){{{C}}_2}(\rho,{\pmb \delta} )}&{{{{A}}_{{k}}}(\rho )} \end{array}} \right], $ 
$ \begin{align} & {{B}_{\text{clp}}}(\rho ,\delta )=\left[\begin{matrix} {{B}_{11}}(\rho ,\delta ) & {{B}_{12}}(\rho ,\delta ) \\ \mathbf{0} & {{B}_{k}}(\rho ) \\ \end{matrix} \right],\\ & {{C}_{\text{clp}}}(\rho ,\delta )=\left[\begin{matrix} {{C}_{11}}(\rho ,\delta ) & \mathbf{0} \\ {{C}_{12}}(\rho ,\delta ) & {{C}_{k}}(\rho ) \\ \end{matrix} \right],\\ & {{D}_{\text{clp}}}(\rho ,\delta )=\mathbf{0},\\ \end{align} $ 
and ${{W}}$ is the matrix associated with the Lyapunov function of the closedloop system,defined here by
$ {{W}} = \left[{\begin{array}{*{20}{c}} {{Y}}&{  ({{Y}}  {{{X}}^{  1}})}\\ {({{Y}}  {{{X}}^{  1}})}&{{{Y}}  {{{X}}^{  1}}} \end{array}} \right]. $ 
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 (LPVC) on nominal plant,SA on 1000 scenarios (whose solution is used as the initial guess of the iterative algorithms), and Wada$'$s algorithm.
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 closedloop 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 closedloop 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 steplike 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 closedloop 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 closedloop system tracks reference model outputs very well. This means that governed by the LPV controller, the system consistently presents a secondorder closedloop characteristic specified by the reference model in a wide speed range (0$\sim$23m/s),which is quite desirable for the highlayer motion planning design. It is further observed from Fig.3 that the variation of vertical speed is kept very small during the highspeed forward flight,which means that the longitudinal phugoid mode is eliminated,demonstrating that the decoupling performance of the controller is satisfactory.
Download:


Fig. 3. Forward and vertical speed responses to steplike excitation on the forward speed channel. 
Download:


Fig. 4. Forward and vertical speed responses to steplike 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 timeconsuming. 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 cuttingplane techniques. For the controller design of a small helicopter$'$s longitudinal motion,a 157dimensional optimization problem,numerical stability would gradually deteriorate,i.e., the condition numbers of the matrices generated from the algorithm become much larger,as the cuttingplanebased 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 largerscale, e.g.,larger than 157dimensions,a more elaborate constraint pruning technique would be necessary.
Ⅴ. CONCLUSIONSIn this paper,we have presented iterative algorithms for the SA. For the robust control of a complex LPV system which is timecostly 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.
[1]  Rugh W J, Shamma J S. Research on gain scheduling. Automatica, 2000, 36(10):14011425 
[2]  Helmersson A. μsynthesis and LFT gain scheduling with real uncertainties. International Journal of Robust and Nonlinear Control, 1988, 8(7):631642 
[3]  Helmerson A. Application of realμ gain scheduling. In:Proceedings of the 35th IEEE Conference on Decision and Control, Kobe, Japan:IEEE, 1996. 16661671 
[4]  Apkarian P, Adams R J. Advanced gainscheduling techniques for uncertain systems. IEEE Transactions on Control System Technology, 1998, 6(1):2132 
[5]  Aouani N, Salhi S, Garcia G, Ksouri M. Robust control analysis and synthesis for LPV systems under affine uncertainty structure. In:Proceedings of the 6th International Multiconference on Systems, Signals and Devices, Djerba, Tunisia:IEEE, 2009. 187191 
[6]  Veenman J, Scherer C W. On robust LPV controller synthesis:a dynamic integral quadratic constraint based approach. In:Proceedings of the 49th IEEE Conference on Decision and Control, Atlanta, USA:IEEE, 2010. 591596 
[7]  Calafiore G, Dabbene F, Tempo R. A survey of randomized algorithms for control synthesis and performance verification. Journal of Complexity, 2007, 23(3):301316 
[8]  Kanev S, De Schutter B, Verhaegen M. An ellipsoid algorithm for probabilistic robust controller design. System & Control Letters, 2003, 49(5):365375 
[9]  Calafiore G, Dabbene F. A probabilistic analytic center cutting plane method for feasibility of uncertain LMIs. Automatica, 2007, 43(12):20222033 
[10]  Calafiore G, Dabbene F, Tempo R. Research on probabilistic methods for control system design. Automatica, 2011, 47(7):12791293 
[11]  Calafiore G, Campi M C. The scenario approach to robust control design. IEEE Transactions on Automatic Control, 2006, 51(5):742753 
[12]  Calafiore G, Fagiano L. Stochastic model predictive control of LPV systems via scenario optimization, Automatica, 2013, 49(6):18611866 
[13]  Campi M C, Garatti S, Prandini M. The scenario approach for systems and control design. Annual Reviews in Control, 2009, 33(2):149157 
[14]  Wang X, Lu G, Zhong Y. Robust H_{∞} attitude control of a laboratory helicopter. Robotics and Autonomous Systems, 2013, 61(12):12471257 
[15]  Liu C, Chen W, Andrews J. Tracking control of smallscale helicopters using explicit nonlinear MPC augmented with disturbance observers. Control Engineering Practice, 2012, 20(3):258268 
[16]  Cai G, Chen B M, Dong X, Lee T H. Design and implementation of a robust and nonlinear flight control system for an unmanned helicopter. Mechatronics, 2011, 21(5):803820 
[17]  Simplício P, Pavel M D, van Kampen E, Chu Q P. An acceleration measurementsbased approach for helicopter nonlinear flight control using incremental nonlinear dynamic inversion. Control Engineering Practice, 2013, 21(8):10651077 
[18]  Wada T, Fujisaki Y. Sequential randomization algorithms:a probabilistic cutting plane technique based on maximum volume ellipsoid center. In:Proceedings of the 2010 IEEE International Symposium on ComputerAided Control System Design, Yokohama, Japan:IEEE, 2010. 15331538 
[19]  Localization methods[Online], available:http://www.stanford.edu/class/ee392o, September 1, 2011 
[20]  Wu F. Control of Linear Parameter Varying Systems[Ph. D. dissertation], University of California at Berkeley, USA, 1995 
[21]  Safonov M, Limbeer D, Chiang R. Simplifying the H_{∞} theory via loop shifting, matrix pencil, and descriptor concepts. International Journal of Control, 1989, 50(6):24672488 
[22]  Calafiore G, Polyak B T. Stochastic algorithms for exact and approximate feasibility of robust LMIs. IEEE Transactions on Automatic Control, 2001, 46(11):17551759 
[23]  Lu B, Wu F. Probabilistic robust linear parametervarying control of an F16 aircraft. AIAA Journal of Guidance, Control, and Dynamics, 2006, 29(6):14541459 
[24]  Heffley R K, Mnich M A. Minimumcomplexity Helicopter Simulation Math Model. Technical Report 19880020435, NASA, USA, 1988 
[25]  Gavrilets V. Autonomous Aerobatic Maneuvering of Miniature Helicopters[Ph. D. dissertation], Massachusetts Institute of Technology, 2003. 
[26]  Löfberg J. YALMIP:A toolbox for modeling and optimization in MATLAB. In:Proceedings of the 2004 IEEE International Symposium on Computer Aided Control Systems Design, Taipei, China:IEEE, 2004. 284289 