control of aircraft lateral axis using mu synthesis -凯发k8网页登录

performance specifications

the illustration below shows a block diagram of the closed-loop system. the diagram includes the nominal aircraft model, the controller k, as well as elements capturing the model uncertainty and performance objectives (see next sections for details).

figure 1: robust control design for aircraft lateral axis

the design goal is to make the airplane respond effectively to the pilot's lateral stick and rudder pedal inputs. the performance specifications include:

hq_p    = 5.0 * tf(2.0,[1 2.0]);
step(hq_p), title('desired response from lateral stick to roll rate (handling quality)')

figure 2: desired response from lateral stick to roll rate.

hq_beta = -2.5 * tf(1.25^2,[1 2.5 1.25^2]);
step(hq_beta), title('desired response from rudder pedal to side-slip angle (handling quality)')

figure 3: desired response from rudder pedal to side-slip angle.

freq = 12.5 * (2*pi);  % 12.5 hz
zeta = 0.5;
yaw_filt = tf(freq^2,[1 2*zeta*freq freq^2]);
lat_filt = tf(freq^2,[1 2*zeta*freq freq^2]);
freq = 4.1 * (2*pi);  % 4.1 hz
zeta = 0.7;
roll_filt = tf(freq^2,[1 2*zeta*freq freq^2]);
aafilters = append(roll_filt,yaw_filt,lat_filt);

from specs to weighting functions

h-infinity design algorithms seek to minimize the largest closed-loop gain across frequency (h-infinity norm). to apply these tools, we must first recast the design specifications as constraints on the closed-loop gains. we use weighting functions to "normalize" the specifications across frequency and to equally weight each requirement.

we can express the design specs in terms of weighting functions as follows:

w_act = ss(diag([1/50,1/20,1/60,1/30]));

w_n = append(0.025,tf(0.0125*[1 1],[1 100]),0.025);
clf, bodemag(w_n(2,2)), title('sensor noise power as a function of frequency')

figure 4: sensor noise power as a function of frequency

w_p = tf([0.05 2.9 105.93 6.17 0.16],[1 9.19 30.80 18.83 3.95]);
clf, bodemag(w_p), title('weight on handling quality spec')

figure 5: weight on handling quality spec.

w_beta = 2*w_p;

here we scaled the weights w_act, w_n, w_p, and w_beta so the closed-loop gain between all external inputs and all weighted outputs is less than 1 at all frequencies.

nominal aircraft model

a pilot can command the lateral-directional response of the aircraft with the lateral stick and rudder pedals. the aircraft has the following characteristics:

the nominal lateral directional model lateralaxis has four states:

these variables are related by the state space equations:

where x = [v; r; p; phi], u = [delta_stab; delta_rud], and y = [beta; p; r; yac].

load lateralaxismodel
lateralaxis

lateralaxis =
 
  a = 
               v         r         p       phi
   v      -0.116    -227.3     43.02     31.63
   r     0.00265    -0.259   -0.1445         0
   p    -0.02114    0.6703    -1.365         0
   phi         0    0.1853         1         0
 
  b = 
        delta_stab   delta_rud
   v        0.0622      0.1013
   r     -0.005252    -0.01121
   p      -0.04666    0.003644
   phi           0           0
 
  c = 
                 v          r          p        phi
   beta     0.2469          0          0          0
   p             0          0       57.3          0
   r             0       57.3          0          0
   yac   -0.002827  -0.007877    0.05106          0
 
  d = 
         delta_stab   delta_rud
   beta           0           0
   p              0           0
   r              0           0
   yac     0.002886    0.002273
 
continuous-time state-space model.

the complete airframe model also includes actuators models a_s and a_r. the actuator outputs are their respective deflection rates and angles. the actuator rates are used to penalize the actuation effort.

a_s = [tf([25 0],[1 25]); tf(25,[1 25])];
a_s.outputname = {'stab_rate','stab_angle'};
a_r = a_s;
a_r.outputname = {'rud_rate','rud_angle'};

accounting for modeling errors

the nominal model only approximates true airplane behavior. to account for unmodeled dynamics, you can introduce a relative term or multiplicative uncertainty w_in*delta_g at the plant input, where the error dynamics delta_g have gain less than 1 across frequencies, and the weighting function w_in reflects the frequency ranges in which the model is more or less accurate. there are typically more modeling errors at high frequencies so w_in is high pass.

% normalized error dynamics
delta_g = ultidyn('delta_g',[2 2],'bound',1.0);
% frequency shaping of error dynamics
w_1 = tf(2.0*[1 4],[1 160]);
w_2 = tf(1.5*[1 20],[1 200]);
w_in = append(w_1,w_2);
bodemag(w_1,'-',w_2,'--')
title('relative error on nominal model as a function of frequency')
legend('stabilizer','rudder','location','northwest');

figure 6: relative error on nominal aircraft model as a function of frequency.

building an uncertain model of the aircraft dynamics

now that we have quantified modeling errors, we can build an uncertain model of the aircraft dynamics corresponding to the dashed box in the figure 7 (same as figure 1):

figure 7: aircraft dynamics.

use the connect function to combine the nominal airframe model lateralaxis, the actuator models a_s and a_r, and the modeling error description w_in*delta_g into a single uncertain model plant_unc mapping [delta_stab; delta_rud] to the actuator and plant outputs:

% actuator model with modeling uncertainty
act_unc = append(a_s,a_r) * (eye(2)   w_in*delta_g);
act_unc.inputname = {'delta_stab','delta_rud'};
% nominal aircraft dynamics
plant_nom = lateralaxis;
plant_nom.inputname = {'stab_angle','rud_angle'};
% connect the two subsystems
inputs = {'delta_stab','delta_rud'};
outputs = [a_s.y ; a_r.y ; plant_nom.y];
plant_unc = connect(plant_nom,act_unc,inputs,outputs);

this produces an uncertain state-space (uss) model plant_unc of the aircraft:

plant_unc

uncertain continuous-time state-space model with 8 outputs, 2 inputs, 8 states.
the model uncertainty consists of the following blocks:
  delta_g: uncertain 2x2 lti, peak gain = 1, 1 occurrences
type "plant_unc.nominalvalue" to see the nominal value and "plant_unc.uncertainty" to interact with the uncertain elements.

analyzing how modeling errors affect open-loop responses

we can analyze the effect of modeling uncertainty by picking random samples of the unmodeled dynamics delta_g and plotting the nominal and perturbed time responses (monte carlo analysis). for example, for the differential stabilizer channel, the uncertainty weight w_1 implies a 5% modeling error at low frequency, increasing to 100% after 93 rad/sec, as confirmed by the bode diagram below.

% pick 10 random samples
plant_unc_sampl = usample(plant_unc,10);
% look at response from differential stabilizer to beta
figure('position',[100,100,560,500]) 
subplot(211), step(plant_unc.nominal(5,1),'r ',plant_unc_sampl(5,1),'b-',10)
legend('nominal','perturbed')
subplot(212), bodemag(plant_unc.nominal(5,1),'r ',plant_unc_sampl(5,1),'b-',{0.001,1e3})
legend('nominal','perturbed')

figure 8: step response and bode diagram.

designing the lateral-axis controller

proceed with designing a controller that robustly achieves the specifications, where robustly means for any perturbed aircraft model consistent with the modeling error bounds w_in.

first we build an open-loop model olic mapping the external input signals to the performance-related outputs as shown below.

figure 9: open-loop model mapping external input signals to performance-related outputs.

to build this model, start with the block diagram of the closed-loop system, remove the controller block k, and use connect to compute the desired model. as before, the connectivity is specified by labeling the inputs and outputs of each block.

figure 10: block diagram for building open-loop model.

% label block i/os
aafilters.u = {'p','r','yac'};    aafilters.y = 'aafilt';
w_n.u = 'noise';                  w_n.y = 'wn';
hq_p.u = 'p_cmd';                 hq_p.y = 'hq_p';
hq_beta.u = 'beta_cmd';           hq_beta.y = 'hq_beta';
w_p.u = 'e_p';                    w_p.y = 'z_p';
w_beta.u = 'e_beta';              w_beta.y = 'z_beta';
w_act.u = [a_s.y ; a_r.y];        w_act.y = 'z_act';
% specify summing junctions
sum1 = sumblk('%meas = aafilt   wn',{'p_meas','r_meas','yac_meas'});
sum2 = sumblk('e_p = hq_p - p');
sum3 = sumblk('e_beta = hq_beta - beta');
% connect everything
olic = connect(plant_unc,aafilters,w_n,hq_p,hq_beta,...
   w_p,w_beta,w_act,sum1,sum2,sum3,...
   {'noise','p_cmd','beta_cmd','delta_stab','delta_rud'},...
   {'z_p','z_beta','z_act','p_cmd','beta_cmd','p_meas','r_meas','yac_meas'});

this produces the uncertain state-space model

olic

uncertain continuous-time state-space model with 11 outputs, 7 inputs, 26 states.
the model uncertainty consists of the following blocks:
  delta_g: uncertain 2x2 lti, peak gain = 1, 1 occurrences
type "olic.nominalvalue" to see the nominal value and "olic.uncertainty" to interact with the uncertain elements.

recall that by construction of the weighting functions, a controller meets the specs whenever the closed-loop gain is less than 1 at all frequencies and for all i/o directions. first design an h-infinity controller that minimizes the closed-loop gain for the nominal aircraft model:

nmeas = 5;		% number of measurements
nctrls = 2;		% number of controls
[kinf,~,gamma_inf] = hinfsyn(olic.nominalvalue,nmeas,nctrls);
gamma_inf

gamma_inf = 0.9700

here hinfsyn computed a controller kinf that keeps the closed-loop gain below 1 so the specs can be met for the nominal aircraft model.

next, perform a mu-synthesis to see if the specs can be met robustly when taking into account the modeling errors (uncertainty delta_g). use the command musyn to perform the synthesis and use musynoptions to set the frequency grid used for mu-analysis.

fmu = logspace(-2,2,60);
opt = musynoptions('frequencygrid',fmu);
[kmu,clperf] = musyn(olic,nmeas,nctrls,opt);

d-k iteration summary:
-----------------------------------------------------------------
                       robust performance               fit order
-----------------------------------------------------------------
  iter         k step       peak mu       d fit             d
    1           5.097        3.487        3.488            12
    2            1.31        1.292        1.312            20
    3           1.242        1.242        1.693            12
    4           1.693        1.544        1.545            16
    5           1.223        1.223        1.551            12
    6           1.533        1.464        1.465            20
    7           1.284        1.283        1.299            12
best achieved robust performance: 1.22

clperf

clperf = 1.2226

here the best controller kmu cannot keep the closed-loop gain below 1 for the specified model uncertainty, indicating that the specs can be nearly but not fully met for the family of aircraft models under consideration.

frequency-domain comparison of controllers

compare the performance and robustness of the h-infinity controller kinf and mu controller kmu. recall that the performance specs are achieved when the closed loop gain is less than 1 for every frequency. use the lft function to close the loop around each controller:

clinf = lft(olic,kinf);
clmu = lft(olic,kmu);

what is the worst-case performance (in terms of closed-loop gain) of each controller for modeling errors bounded by w_in? the wcgain command helps you answer this difficult question directly without need for extensive gridding and simulation.

% compute worst-case gain as a function of frequency
opt = wcoptions('varyfrequency','on');
% compute worst-case gain (as a function of frequency) for kinf
[mginf,wcuinf,infoinf] = wcgain(clinf,opt);
% compute worst-case gain for kmu
[mgmu,wcumu,infomu] = wcgain(clmu,opt);

you can now compare the nominal and worst-case performance for each controller:

clf
subplot(211)
f = infoinf.frequency;
gnom = sigma(clinf.nominalvalue,f);
semilogx(f,gnom(1,:),'r',f,infoinf.bounds(:,2),'b');
title('performance analysis for kinf')
xlabel('frequency (rad/sec)')
ylabel('closed-loop gain');
xlim([1e-2 1e2])
legend('nominal plant','worst-case','location','northwest');
subplot(212)
f = infomu.frequency;
gnom = sigma(clmu.nominalvalue,f);
semilogx(f,gnom(1,:),'r',f,infomu.bounds(:,2),'b');
title('performance analysis for kmu')
xlabel('frequency (rad/sec)')
ylabel('closed-loop gain');
xlim([1e-2 1e2])
legend('nominal plant','worst-case','location','southwest');

the first plot shows that while the h-infinity controller kinf meets the performance specs for the nominal plant model, its performance can sharply deteriorate (peak gain near 15) for some perturbed model within our modeling error bounds.

in contrast, the mu controller kmu has slightly worse performance for the nominal plant when compared to kinf, but it maintains this performance consistently for all perturbed models (worst-case gain near 1.25). the mu controller is therefore more robust to modeling errors.

time-domain validation of the robust controller

to further test the robustness of the mu controller kmu in the time domain, you can compare the time responses of the nominal and worst-case closed-loop models with the ideal "handling quality" response. to do this, first construct the "true" closed-loop model clsim where all weighting functions and hq reference models have been removed:

kmu.u = {'p_cmd','beta_cmd','p_meas','r_meas','yac_meas'};
kmu.y = {'delta_stab','delta_rud'};
aafilters.y = {'p_meas','r_meas','yac_meas'};
clsim = connect(plant_unc(5:end,:),aafilters,kmu,{'p_cmd','beta_cmd'},{'p','beta'});

next, create the test signals u_stick and u_pedal shown below

time = 0:0.02:15;
u_stick = (time>=9 & time<12);
u_pedal = (time>=1 & time<4) - (time>=4 & time<7);
clf
subplot(211), plot(time,u_stick), axis([0 14 -2 2]), title('lateral stick command')
subplot(212), plot(time,u_pedal), axis([0 14 -2 2]), title('rudder pedal command')

you can now compute and plot the ideal, nominal, and worst-case responses to the test commands u_stick and u_pedal.

% ideal behavior
idealresp = append(hq_p,hq_beta);
idealresp.y = {'p','beta'};
% worst-case response
wcresp = usubs(clsim,wcumu);
% compare responses
clf
lsim(idealresp,'g',clsim.nominalvalue,'r',wcresp,'b:',[u_stick ; u_pedal],time)
legend('ideal','nominal','perturbed','location','southeast');
title('closed-loop responses with mu controller kmu')

the closed-loop response is nearly identical for the nominal and worst-case closed-loop systems. note that the roll-rate response of the aircraft tracks the roll-rate command well initially and then departs from this command. this is due to a right-half plane zero in the aircraft model at 0.024 rad/sec.