www.gusucode.com > robust 案例源码程序 matlab代码 > robust/MIMORobustnessAnalysisExample.m

    %% MIMO Robustness Analysis
% You can create and analyze uncertain state-space models made up of uncertain
% state-space matrices. In this example, create a MIMO system with
% parametric uncertainty and analyze it for robust stability and worst-case
% performance.
% 
%%
% Consider a two-input, two-output, two-state system whose
% model has parametric uncertainty in the state-space matrices. First
% create an uncertain parameter |p|.  Using the parameter, make uncertain
% |A| and |C| matrices.  The |B| matrix happens to be not-uncertain,
% although you will add frequency-domain input uncertainty to the model
% later. 
p = ureal('p',10,'Percentage',10);
A = [0 p;-p 0];
B = eye(2);
C = [1 p;-p 1];
H = ss(A,B,C,[0 0;0 0])
%%
% You can view the properties of the uncertain system |H| using the |get|
% command. 
get(H)
%%
% Most properties behave in the same way as the corresponding properties of
% |ss| objects.  The property |NominalValue| is itself an |ss| object.
%% Adding Independent Input Uncertainty to Each Channel
% The model for |H| does not include actuator dynamics. Said differently,
% the actuator models are unity-gain for all frequencies. 
%
% Nevertheless, the behavior of the actuator for channel 1 is modestly
% uncertain (say 10%) at low frequencies, and the high-frequency behavior
% beyond 20 rad/s is not accurately modeled. Similar statements hold for
% the actuator in channel 2, with larger modest uncertainty at low
% frequency (say 20%) but accuracy out to 45 rad/s.
%
% Use |ultidyn| objects |Delta1| and |Delta2| along with shaping filters
% |W1| and |W2| to add this form of frequency domain uncertainty to the
% model. 
W1 = makeweight(.1,20,50);
W2 = makeweight(.2,45,50);
Delta1 = ultidyn('Delta1',[1 1]);
Delta2 = ultidyn('Delta2',[1 1]);
G = H*blkdiag(1+W1*Delta1,1+W2*Delta2)
%%
% Note that |G| is a two-input, two-output uncertain system, with
% dependence on three uncertain elements, |Delta1|, |Delta2|, and |p|.  It
% has four states, two from |H| and one each from the shaping filters |W1|
% and |W2|, which are embedded in |G|.  
%
% You can plot a 2-second step response of several samples of |G| The 10%
% uncertainty in the natural frequency is obvious. 
stepplot(G,2)
%%
% You can create a Bode plot of samples of |G|. The high-frequency
% uncertainty in the model is also obvious. For clarity, start the Bode
% plot beyond the resonance.
bodeplot(G,{13 100})
%% Closed-Loop Robustness Analysis
% Load the controller and verify that it is two-input and two-output.
load(fullfile(matlabroot,'examples','robust','mimoKexample.mat'))
size(K)
%%
% You can use the command |loopsens| to form all the standard
% plant/controller feedback configurations, including sensitivity and
% complementary sensitivity at both the input and output. Because |G| is
% uncertain, all the closed-loop systems are uncertain as well.
F = loopsens(G,K)
%%
% |F|  is a structure with many fields. The poles of the nominal
% closed-loop system are in |F.Poles|, and |F.Stable| is 1 if the nominal
% closed-loop system is stable. In the remaining 10 fields, |S| stands for
% sensitivity, |T| or complementary sensitivity, and |L| for open-loop
% gain. The suffixes |i| and |o| refer to the input and output of the
% plant.  Finally, |P| and |C| refer to the plant and controller.
%
% Hence, |Ti| is mathematically the same as: 
%%
% 
% $$K(I+GK)^{-1}G$$
% 
%% 
% |Lo| is |G*K|, and |CSo| is mathematically the same as
%%
% 
% $$K(I+GK)^{-1}$$
% 
%% 
% Examine the transmission of disturbances at the plant input to the plant
% output by plotting responses of |F.PSi|. Graph some samples along with
% the nominal.
bodemag(F.PSi.NominalValue,'r+',F.PSi,'b-',{1e-1 100})

%% Nominal Stability Margins
% You can use |loopmargin| to investigate loop-at-a-time gain and phase
% margins, loop-at-a-time disk margins, and simultaneous multivariable
% margins. They are computed for the nominal system and do not reflect the
% uncertainty models within |G|.
%
% Explore the simultaneous margins individually at the plant input,
% individually at the plant output, and simultaneously at both input and output. 
[I,DI,SimI,O,DO,SimO,Sim] = loopmargin(G,K);
%%
% The third output argument is the simultaneous gain and phase variations
% allowed in all input channels to the plant.
SimI
%%
% This information implies that the gain at the plant input can vary in
% both channels independently by factors between (approximately) 1/8 and 8,
% as well as phase variations up to 76 degrees.
%
% The sixth output argument is the simultaneous gain and phase variations
% allowed in all output channels to the plant.
SimO
%% 
% Note that the simultaneous margins at the plant output are similar to
% those at the input. This is not always the case in multiloop feedback
% systems.
%
% The last output argument is the simultaneous gain and phase variations
% allowed in all input and output channels to the plant. As expected, when
% you consider all such variations simultaneously, the margins are somewhat
% smaller than those at the input or output alone.
%
Sim
%%
% Nevertheless, these numbers indicate a generally robust closed-loop
% system, able to tolerate significant gain (more than +/-50% in each
% channel) and 30 degree phase variations simultaneously in all input and
% output channels of the plant.

%% Robust Stability Margin
% With |loopmargin|, you determined various stability margins of the nominal
% multiloop system. These margins are computed only for the nominal system
% and do not reflect the uncertainty explicitly modeled by the |ureal| and
% |ultidyn| objects.  When you work with a detailed uncertainty
% model, the conventional stability margins computed by |loopmargin| 
% may not accurately reflect how close the system is from being unstable.
% You can then use |robstab| to compute the robust stability margin for
% the specified uncertainty.
% 
% In this example, use |robstab| to compute the robust stability margin for
% the uncertain feedback loop comprised of |G| and |K|. You can use
% any of the closed-loop transfer functions in |F = loopsens(G,K)|.  All of
% them, |F.Si, F.To|, etc., have the same internal dynamics, and hence
% their stability properties are the same.  
opt = robOptions('Display','on');
stabmarg = robstab(F.So,opt)
%%
% This analysis confirms what the |loopmargin| analysis suggested. The
% closed-loop system is quite robust, in terms of stability,  to the
% variations modeled by the uncertain parameters |Delta1|, |Delta2|, and
% |p|.  In fact, the system can tolerate more than twice the modeled
% uncertainty without losing closed-loop stability.

%% Worst-Case Gain Analysis
% You can plot the Bode magnitude of the nominal output sensitivity
% function. It clearly shows decent disturbance rejection in all channels
% at low frequency.
bodemag(F.So.NominalValue,{1e-1 100})
%%
% You can compute the peak value of the maximum singular value of the
% frequency response matrix using |norm|.  
[PeakNom,freq] = getPeakGain(F.So.NominalValue)
%% 
% The peak is about 1.13. What is the
% maximum output sensitivity gain that is achieved when the uncertain
% elements |Delta1|, |Delta2|, and |p| vary over their ranges? You can use
% |wcgain| to answer this.
[maxgain,wcu] = wcgain(F.So);
maxgain
%% 
% The analysis indicates that the worst-case gain is somewhere between 2.1
% and 2.2. The frequency where the peak is achieved is about 8.5.
% 
% Use |usubs| to replace the values for |Delta1|, |Delta2|, and |p| that achieve
% the gain of 2.1. Make the substitution in the output
% complementary sensitivity, and do a step response.
step(F.To.NominalValue,usubs(F.To,wcu),5)
%%
% The perturbed response, which is the worst combination of uncertain
% values in terms of output sensitivity amplification, does not show
% significant degradation of the command response. The settling time is
% increased by about 50%, from 2 to 4, and the off-diagonal coupling is
% increased by about a factor of about 2, but is still quite small.
%%
% You can also examine the worst-case frequency response alongside the
% nominal and sampled systems using |wcsigma|.
wcsigma(F.To,{1e-1,100})