_ Exercice 1 SISO systeem - DOC by decree

VIEWS: 48 PAGES: 13

									Exercise 18 WB 2421 2004: SISO system
IMPORTANT: To understand how a particular Matlab command works, type help followed by the
name of the command. Example: help plot. USE M.file hinf_siso.m

PART 1: Representation of systems in Matlab.

One way to represent a Single Input Single Output system is through its transfer function, i.e. given its
numerator and denominator polynomials. A polynomial in Matlab is represented by the vector of its
coefficients, and the multiplication of polynomials is done with the command conv.

Example: To represent the system

                               s 2  10s  13
                   P( s ) 
                              ( s  1 )( s  2 )

use the commands:

num=[1 10 13];                             % numerator of P
den=conv([1 1],[1 2]);                     % denominator of P
P=nd2sys(num,den);

The command nd2sys converts a SISO transfer function representation into the system matrix, which is
the structure used in the Mu-Toolbox to represent a system.

Now you do: Consider the model of the motion system:
                                                                                                        x2


                                                                                         c
                                                                          F
                                                                                 M1                M2


    H
         x2
             2
                               b g
                         ds  c  0.02                                                   d

                 b     2
                              g b g b
         F s J1  J 2 s  J1  J 2 ds  J1  J 2 c             g
where J1  J 2  4.8e  6         c  0.22         d  1.e  4 .

Build the system G in Matlab, determine the zero-pole-gain form, the transfer function form, and the
state-space realization of G.
numG=0.02*([d c]);
denG=[(J1*J2) d*(J1+J2) c*(J1+J2) 0 0];
Gmu=nd2sys(numG,denG);                                             % representation according to the Mu-
                                                                   % Tools

G=tf(numG,denG)                                                    % transfer function form of G
ZP_G=zpk(tf_G)                                                     % zero-pole-gain form of G



                                                              1
[A,B,C,D]=unpck(G)                           % state-space realization of G

% Calculate the poles and the zeros and examine the properties of the system G

plant_poles=spoles(G)                   % poles of the system
plant_zeros=szeros(G)                   % zeros of the system

minfo(G)                                % Returns structure information about G

seesys(G)                               % Returns the state-space realization of G

PART 2: Calculate and plot frequency responses of the system G.
In order to calculate a frequency response of a system, first a frequency vector
can be defined. This is done with the command LOGSPACE. Once the frequency vector
has been defined, the frequency response can be calculated with the command FRSP.
Notice that this command requires a frequency vector in rad/sec and not in Hertz!
Also the commands bode and nyquist can be used directly!

% Perform the following:

hz=logspace(-5,5,200);                            %   Frequency vector in Hz of 200
                                                                    5       5
                                                  %   points from 10 Hz to 10 Hz
wz=2*pi*hz;                                       %   Frequency vector in rad/s
ghat=frsp(G,wz);                                  %   varying frequency response
                                                  %   matrix of G

figure(1);
vplot('bode',ghat);grid on;                       % Bode plot of G
subplot(211);grid on;title('Bode plot of G');

% Now suppose that we close the loop around G with a unitary negative feedback. Can
 you say whether the closed loop system is stable or not?


                     ref   +                                     out
                                              G
                               -




To do this plot de Nyquist plot of the open loop system G and check the Nyquist
stability requirement

figure(2)
subplot(211);
vplot('nyq',ghat);grid on;                % Nyquist plot of G
title('Nyquist plot of G');



subplot(212);


                                         2
nyquist(tf_G,{1,1000});grid on              %   plots the Nyquist plot between
                                            %   1 [rad/s] and 1000 [rad/s].
                                            %   zoom on to see a part of the plot near
                                            %   the (–1)-point.


PART 3: Calculate and plot time response of the system G
In order to calculate a time response of a system, first a time vector and an input
signal should be defined. The time can be defined with the command LINSPACE. The
input signal is defined with the command SIGGEN. Then the time response of the
system for the defined input is obtained using the command TRSP and can be plotted
using VPLOT with the option 'iv,m'.
Also the commands step, gensig and lsim can be used.

% Plot the response of G to a step input:

t=linspace(0,5,2000);                            %   Time vector of 2000 points from
                                                 %   0 to 5 sec.
ustep=siggen('1',t);                             %   Define step input.
ystep=trsp(G,ustep,5);                           %   step response of G within 5[sec]

figure(3)
subplot(211);
vplot('iv,m',ystep);grid on;
title('step response')

% Plot the response of G to a sine of frequency of 5 Hz:

usin=siggen('sin(2*pi*5*t)',t);                      %Define sine input.
ysin=trsp(G,usin,5);                                 %sine response of G

subplot(212);
vplot('iv,m',ysin);grid on;
title('response to sinusoidal input')

PART 4: Analysis of the closed loop system
Consider the classical closed loop scheme of figure 1.



             ref    +                                                  out
                                    K                         G
                          -




                         Figure 1

Now you do: Based on the analysis of the Bode and Nyquist plots of PART 2, try to
determine a PD controller K for which the closed loop system is stable with a phase




                                            3
margin of at least 30 [deg] and with a peak of the sensitivity function not bigger
than 2.
To do this, first determine the interconnected system in open loop form of figure
2, and then close the loop with the controller.
                           z                                z2
                              1




              w      +            y        u
                                                     G
                         _




                             z1                            w
                             z2             P
                             y                             u


                             Figure 2: generalized plant structure

% Import the following to build the interconnection:


systemnames='Gmu';                    % Defines the system blocks in the structure
                                      % of figure 2

inputvar='[w;u]';                     % Defines the input variables in the structure

outputvar='[w-Gmu;Gmu;w-Gmu]';        % Defines the output variables z1,z2,y in the
                                      % structure

input_to_Gmu='[u]';                   % Specifies the input of the plant G
sysoutname='Pmu';                     % Gives the name of the interconnection

cleanupsysic='yes';
sysic;                                % Runs the interconnection routine

seesys(Pmu)                           % Returns the structure of Pmu
minfo(Pmu)

                                    % Convert to lti format
[a,b,c,d]=unpck(Pmu);
P=ss(a,b,c,d); P.InputName={'w';'u'}; P.OutputName={'z1';'z2';'y'};




                                           4
% Close the loop with a PD controller K as shown in figure 4.


                            z1                                  w
                            z2                   P
                            y                                    u




                                                 K

                                 %Figure 4: closed loop system.
% Import the following:

                                   % Put here your choice of K=Kp+Kv*s
%Kp=?;
%Kv=?;

Kp=0.004;
Kv=0.00016;
numk=[Kv Kp];
denk=[1.e-8 1];
K=tf(numk,denk);
Pcl=lft(P,K);                           % Closed loop system matrix

                                        %   the commando lft is the Linear Fractional
                                        %   transformation of the plant P with controller K
                                        %   Type "help lft" for more information.
                                        %   For our specific interconnection we have
                                        %   where S is sensitivity and T is
                                        %   complementary sensitivity

% Calculate the poles of the closed loop system:

         % Is the closed loop stable?

 closed_loop_poles=eig(Pcl)


% Determine and plot the frequency responses of the following functions of the
closed loop system:

          The complementary sensitivity function: T
          The sensitivity function: S

T=Pcl('z2','w');                   % The commando selects subsystems from a Pcl




                                                 5
S=Pcl('z1','w');


That=frsp(T,wz);                                  % frequency response of T
Shat=frsp(S,wz);                                  % frequency response of S

figure(4);
subplot(311);
vplot('liv,lm',That);grid on;
title('Complementary sensitivity: T')

subplot(312);grid on;
vplot('liv,lm',Shat);grid on;
title('Sensitivity: S')

% Determine the peak values of the sensitivity and the complementary sensitivity
functions:

Tpeak=norm(T,inf)
Speak=norm(S,inf)

% Determine and plot the step response of the closed loop system. Is there a
steady-state error?

yclstep=trsp(T,ustep,5);                                 % step response of the closed loop

subplot(313);
vplot('iv,m',yclstep);grid on;
title('closed_loop step response')


PART 5: Robust stability analysis for the PD controller.
Remark: The H-infinity norm is a suitable measure to assess robustness
        of a control system. In Matlab it can be calculated using the
        command HINFNORM or NORM.

Suppose that the plant is not perfectly known and that we can model it as the
nominal value G with an output multiplicative uncertainty as:

                                  b               g
                     G p ( s )  1  Wu ( s )( s ) G( s )

where Wu ( s ) gives the profile of the uncertainty as a function of the frequency.

Suppose:

                     Wu ( s ) 
                                     b
                                  10 s  0.01 g
                                    s  100

Determine whether the controlled system is robustly stable against the defined
uncertainty?

For robust stability the following H-infinity condition should be satisfied:




                                                     6
                     Wu T   
                                 1


The following block scheme gives the weighted plant without controller:

                                                                    z2


                            z1
                                                                  Wu

        w        +                    y       u
                                                        G
                     _



                     Figure 5



                                          ~
% Determine the weighted plant P of figure 5:

% Import the weight Wu

numwu=10*[1 0.01];
denwu=[1 100];
wu=nd2sys(numwu,denwu);



% Determine the weighted complementary sensitivity function and its H  -norm:


norm(wu*T,inf)




figure(5)

subplot(211);
vplot('liv,lm',That);grid on;hold on;
title('Unweighted complementary sensitivity & Inverse weighting')
vplot('liv,lm',wu_invhat,'r');grid on;

subplot(212);
vplot('liv,lm',wuThat);grid on;
title('Weighed complementary sensitivity wuT')

                                                  % Is the controlled system robustly
                                                  % stable?



                                              7
                                                         % Are the plots in agreement with the
                                                         % H-infinity norm?



PART 6: Robust performance analysis for the PD controller.
Consider the uncertain system of PART5 and assume that the performance of the
closed loop system is specified in terms of the weighted sensitivity function. Take
as the performance weight Wp the inverse of the sensitivity function computed in
PART4.    Is the closed loop robustly performing for every uncertainty such that
    1?

For robust performance the following H-infinity condition should hold:

                     Wu T
                                  1.
                     Wp S
                            


The following block scheme gives the weighed plant with Wp and Wu:

                            z1                                               z2



                         Wp                                                 Wu

        w      +                    y             u
                                                               G
                     _



                    Figure 6




% Determine the weighted plant tP of figure 6, close the loop with the constant
% controller, and Determine the H  norm of the controlled sytem:

wp=minv(S);
wp_inv=S;
wp_invhat=Shat;
Wm=daug(wp,wu,1);                       %   The   daug   command   groups   systems   in   diagonal
form.
tP=mmult(Wm,P);
seesys(tP)

% Close the loop




                                                  8
tN=starp(tP,K);
seesys(tN);

tNhat=frsp(tN,wz);

hinfnorm(tNhat)                            % Is the closed loop system robustly performing?



PART 7: S-design
The performance that can be achieved by a simple constant controller is very poor.
For instance the bandwidth of the closed loop system is very low. Shaping the
sensitivity function S(s) by H  control is a suitable design tool to improve
performance. Typical specifications that are related to the sensitivity function
are the bandwidth of the closed loop system, the tracking error and disturbance
attenuation.

Consider a performance weighting function of the form:


                                   2                 2
                                  s  2ns  n
                           Wp 
                                    ( s  a ) ( s  b)

Select appropriate parameters  n ,  , a                    and b in order to satisfy the following
specifications:

1. Bandwidth between 5 and 15 rad/s. (as an example)
2. Disturbance attenuation at low frequencies.
3. Zero steady-state error of the step response.

Perform an S-design that achieves the above specifications:
The block scheme on which the S-design is based is given in figure 7:

                            z1

                                                                                z2
                           Wp

        w         +                    y                 u
                                                                      G
                       _




                      Figure 7

Import the following in Matlab:

%select the appropriate values of the parameters:



                                                         9
wb=?;                                           % Desired bandwidth
z=?;
a=?;
b=?;

numwp=[1 2*z*wb wb^2];
denwp=conv([1 a],[1 b]);
Wp=tf(numwp,denwp);

% Plot the weight function Wp:

wphat=frsp(wp,w);

figure(6);

vplot('liv,lm',wphat);grid on;
title('Weight function Wp')

% Determine the weighted plant tP

W=append(Wp,1,1);

WP= W*P;
[a,b,c,d]=ssdata(WP); WPmu=pck(a,b,c,d);
seesys(WPmu)

seesys(tP)

% Given the weighted plant WP, compute the Hinf controller K(s) which stabilizes
the plan tP and which yields a closed-loop peak smaller than gopt. gopt is the
smallest achievable peak value:

[gopt,Kmu] = hinflmi(MPmu,[1 1],2,1e-3);

                                       % The command hinflmi computes the
                                       % H-infinity contoller according to the
                                       % LMI-Tools.

                                       %   tP   : Weighted open-loop system
                                       %   [1,1]: 1=Number of measured variables
                                       %          1=Number of controlled signals
                                       %   2    : Try to achieve peak value smaller
                                       %          than 2
                                       %   1e-3 : Accuracy

% Determine the state-space realization of the controller K, the transfer function,
% the zero-pole-gain form, and the poles of K.

[Ak,Bk,Ck,Dk]=unpck(Kmu);                          % state-space realization
K=ss(Ak,Bk,Ck,Dk);
[numk,denk]=SS2TF(Ak,Bk,Ck,Dk);                  % numerator and denominator
tf_k=tf(numk,denk)
ZP_kS=zpk(tf_k)
K_poles=spoles(K)

% Plot the magnitude of the controller and interpret what kind of action it


                                           10
% performs. Plot also the magnitudes of the model G, and of the open loop system
GK:

Khat=frsp(K,wz);
GK=mmult(G,K);
GKhat=frsp(GK,wz);

figure(7)

vplot('liv,lm',ghat); grid on;hold on;
vplot('liv,lm',Khat,'r');hold on
vplot('liv,lm',GKhat,'b');
title('Model G, controller K, open loop GK : S-design')


% Analysis of controlled system

% Close the loop of the plant P (not weighted) with the H-infinity controller K.
% Calculate the closed loop poles

Pcl2=lft(P,K);
seesys(Pcl2)
cl_poles=spoles(Pcl2)

% Plot the frequency response of the following:

%    *   Sensitivity function S
%    *   Inverse of weight Wp
%    *   Complementary sensitivity T
%    *   Inverse of weight Wu
%    *   Weighted sensitivity WpS
%    *   Weighted complementary sensitivity
%    *   Achieved performance gopt




% Is the controlled system robustly stable against the multiplicative
% uncertainty of PART 5?

norm(Wu*Pcl2('z2','w'),inf)

% Plot the step response of the closed loop system:

figure(9);grid on
step(Pcl('z2','w'),5)
title('closed_loop step response: S-design')          %   Is the steady state error
                                                      %   equal to zero?
                                                      %   Is the steady state error
                                                      %   equal to zero?
PART 8: S/T-design:
The mixed S/T design allows to combine performance requirements, expresse d in terms
of the sensitivity function, and robustness requirements, expressed in terms of the
complementary sensitivity function. In this design, the objective is the
minimization of:



                                              11
                          F SI
                          W
                          GT J
                          HK
                          W
                               p

                               u       


where Wu reflects the size of the multiplicative uncertainty. Assume Wu as given in
PART 5 and Wu as in PART 7.

Perform an S/T-design that achieves the specifications in PART 7:
The mixed S/T-design is based on the following scheme:


                          z1                                       z2



                          Wp                                     Wu

        w      +                   y       u
                                                       G
                      _




                   Figure(8)

% Determine the weighted plant WmPmu with Wp and Wu as in figure 8

Wm=append(Wp,Wu,1);
WmP=Wm*P;
[a,b,c,d]=ssdata(WmP); WmPmu=pck(a,b,c,d);
seesys(WmPmu)

% Compute the Hinfinity controller K

[gopt,K2mu] = hinflmi(WmPmu,[1 1],2,1e-3);

[Ak2,Bk2,Ck2,Dk2]=unpck(K2mu);                  % state-space realization
K2=ss(Ak2,Bk2,Ck2,Dk2);

[numk2,denk2]=tfdata(tf(K2));                   % numerator and denominator
[K2_poles,K2_zeros]=pzmap(K2)



% Plot the magnitude of the controller and interpret what kind of action it
 performs. Plot also the magnitudes of the model G, and of the open loop system GK:

Khat=frsp(K2,wz);
GK2=mmult(G,K2);
GKhat=frsp(GK2,wz);



                                           12
figure(10)

vplot('liv,lm',ghat); grid on;hold on;
vplot('liv,lm',Khat,'r');hold on
vplot('liv,lm',GKhat,'b');
title('Model G, controller K, open loop Gk : S/T-design')


       % What is the difference between the controller of the S-design and of
       % the S/T-design? ( use e.g. magnitude plots and/or the zero-pole form of
       % the controllers)


%   Analysis of controlled system
% Close the loop of the plant P (not weighted) with the H-infinity controller K.
% Calculate the closed loop poles

Pcl3=lft(P,K2);
seesys(Pcl3)
cl_poles=spoles(Pcl3)


% Plot the frequency response of the following:
%    * Sensitivity function S
%    * Inverse of weight Wp
%    * Complementary sensitivity T
%    * Inverse of weight Wu
%    * The weighted sensitivity WpS
%    * The weighted complementary sensitivity WuT
%    * The achieved performance gopt




% Plot the step response of the closed loop system:



% Examine robust stability and robust performance

% Compare the performance and robustness aspects of S-design and S/T-design


% Plot the frequency response of KS, the transfer function from the reference
% signal to the output of the controller.


% Draw conclusions




                                         13

								
To top