%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% file uebung_control2_e.m
%
% Exercise to lesson MATLAB Simulink - Control System TB II
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

echo on;
format compact;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creation of LTI-models
uebung_control1_e

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Analysing LTI-models
% Model dynamics 

class(pt1)              % Type of LTI-model                                  
dcgain(pt1)             % DC gain value                                      
damp(pt1)               % natural frequencies and damping                    
[p,z] = pzmap(pt1)      % Poles and zeros                                    
allmargin(pt1)          % Gain margin, Phasen margin, delay margin, stability

pause    
class(int)              % Type of LTI-model                                  
dcgain(int)             % DC gain value                                      
damp(int)               % natural frequencies and damping                    
[p,z] = pzmap(int)      % Poles and zeros                                    
allmargin(int)          % Gain margin, Phasen margin, delay margin, stability

pause
class(reg)              % Type of LTI-model                                  
dcgain(reg)             % DC gain value                                      
damp(reg)               % natural frequencies and damping                    
[p,z] = pzmap(reg)      % Poles and zeros                                    
allmargin(reg)          % Gain margin, Phasen margin, delay margin, stability

pause
class(sysfrd)           % Type of LTI-model                                  
allmargin(sysfrd)       % Gain margin, Phasen margin, delay margin, stability
                           
pause                      
                           
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time responses 

figure , set(gcf,'Name','pt1, dpt1')
subplot(221), impulse(pt1,dpt1)       % Impulse response
subplot(222), step(pt1,dpt1)          % Step response 
subplot(223), pzmap(pt1)              % Zero-Pole-Map
subplot(224), pzmap(dpt1)             % Zero-Pole-Map

pause
figure , set(gcf,'Name','int, dint')
subplot(221), impulse(int,dint)       % Impulse response               
subplot(222), step(int,dint)          % Step response                  
subplot(223), pzmap(int)              % Zero-Pole-Map                  
subplot(224), pzmap(dint)             % Zero-Pole-Map                  

pause
figure , set(gcf,'Name','reg, dreg')
subplot(221), impulse(reg,dreg)       % Impulse response               
subplot(222), step(reg,dreg)          % Step response                  
subplot(245), pzmap(reg)              % Zero-Pole-Map                  
subplot(246), pzmap(dreg)             % Zero-Pole-Map                  
subplot(224), initial(reg,dreg,1.4)   % Initial response 

pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Frequency responses

figure , set(gcf,'Name','pt1, dpt1')
subplot(221), bode(pt1,dpt1)          % Bode diagram
subplot(222), margin(pt1)             % Bode margin plot

pause
figure , set(gcf,'Name','int, dint')
subplot(221), bode(int,dint)          % Bode diagram                 
subplot(222), margin(int)             % Bode margin plot             

pause
figure , set(gcf,'Name','reg, dreg')
subplot(221), bode(reg,dreg)          % Bode diagram                 
subplot(222), margin(reg)             % Bode margin plot             

pause
figure , set(gcf,'Name','sysfrd')
subplot(221), bode(sysfrd)            % Bode diagram                 
subplot(222), margin(sysfrd)          % Bode margin plot             

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compensator systems
% Analyzing model properties

g = - forward * mess

pause
dcgain(plant)                         % DC gain value                                      
damp(plant)                           % natural frequencies and damping                    
[p,z] = pzmap(plant)                  % Poles and zeros                                    
allmargin(plant)                      % Gain margin, Phasen margin, delay margin, stability
figure , set(gcf,'Name','plant')
subplot(221), step(plant)             % Step response   
subplot(222), pzmap(plant)            % Zero-Pole-Map   
subplot(223), margin(plant)           % Bode margin plot
%subplot(224), nyquist(plant)          % Nyquist plot

pause
dcgain(forward)                       % DC gain value                                      
damp(forward)                         % natural frequencies and damping                    
[p,z] = pzmap(forward)                % Poles and zeros                                    
allmargin(forward)                    % Gain margin, Phasen margin, delay margin, stability
figure , set(gcf,'Name','forward')                                                         
subplot(221), step(forward)           % Step response                                      
subplot(222), pzmap(forward)          % Zero-Pole-Map                                      
subplot(223), margin(forward)         % Bode margin plot                                   
%subplot(224), nyquist(forward)        % Nyquist plot

pause
dcgain(g)                             % DC gain value                                      
damp(g)                               % natural frequencies and damping                    
[p,z] = pzmap(g)                      % Poles and zeros                                    
allmargin(g)                          % Gain margin, Phasen margin, delay margin, stability
figure , set(gcf,'Name','g')                                                               
subplot(221), step(g)                 % Step response                                      
subplot(222), pzmap(g)                % Zero-Pole-Map                                      
subplot(223), margin(g)               % Bode margin plot                                   
%subplot(224), nyquist(g)              % Nyquist plot

pause
dcgain(sys)                           % DC gain value                                      
damp(sys)                             % natural frequencies and damping                    
[p,z] = pzmap(sys)                    % Poles and zeros                                    
allmargin(sys)                        % Gain margin, Phasen margin, delay margin, stability
dcgain(sys2)                          % DC gain value                                      
damp(sys2)                            % natural frequencies and damping                    
[p,z] = pzmap(sys2)                   % Poles and zeros                                    
allmargin(sys2)                       % Gain margin, Phasen margin, delay margin, stability
figure , set(gcf,'Name','sys, sys2')
subplot(221), step(sys,sys2,'r')      % Step response                                      
subplot(222), pzmap(sys,sys2,'r')     % Zero-Pole-Map                                      
subplot(223), bode(sys,sys2,'r')      % Bode margin plot                                   
%subplot(224), nyquist(sys,sys2,'r')   % Nyquist plot


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Different input signals 

[u,t] = gensig('square',2) ;
figure
subplot(221), lsim(sys,u,t)
subplot(222), lsim(sys2,u,t)

[u,t] = gensig('sin',1) ;
subplot(223), lsim(sys,u,t)
subplot(224), lsim(sys2,u,t)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Variation of compensator time constant
%% a) Create LTI-models
sys  = tf(1,[.1 1])                         % LTI-model sys
% or
s = tf('s')
sys  = 1 / (1+0.1*s)                        % alternative LTI-model sys

%% b) Create P-Compensator 
regP   = tf(2,1)                            % P-Compensator 

%% c) Transfer function of closed feedback loop
regsysP   = feedback(sys*regP,1)            

%% d) Step responses
figure
step(sys,'-',regsysP,':')
legend('sys','regsysP')

%% e) Create PI-Compensator 
T_R     = 0.01 ;
V_R     = 2 ;
regPI   = tf(V_R*[T_R 1],[T_R 0])             % PI-Compensator 

%% f) Transfer function of closed feedback loop
regsysPI  = feedback(sys*regPI,1)           

%% g) Step response PI-compensated plant:
%axis([0 3.5 0 1])
step(sys,'b-',regsysP,'b:',regsysPI,'b-.')
legend('sys','regsysP','regsysPI')

%% h) Variation compensator time constant T_R
% Appropriate solution
step(sys,'k:')
T_R = [0.01 0.1 0.2 ]
for n=1:length(T_R)
    regPI    = tf(V_R*[T_R(n) 1],[T_R(n) 0]);         % Create PI-Compensator 
    regsysPI = feedback(sys*regPI,1)           
    damp(regsysPI);
    [p,z] = pzmap(regsysPI)
    hold on
    step(regsysPI)   
    pause
end

% Solution for sophisticated output 
step(sys,'r-')
lt = {} ;
T_R = [0.01 0.1 0.2 ]
for n=1:length(T_R)
    regPI    = tf(2*[T_R(n) 1],[T_R(n) 0]);         % Create PI-Compensator 
    regsysPI = feedback(sys*regPI,1)           
    damp(regsysPI);
    [p,z] = pzmap(regsysPI)
    hold on
    stil = {'b--' 'k:' 'g-.'} ;                     
    step(regsysPI,stil{1+mod(n-1,length(T_R))})     
    lt = { lt{:} , ['sysregPI - T_R(',num2str(n),')=',num2str(T_R(n))] }';
    legend('sys' , lt{:} , 'Location','SouthEast')
    axis([0 0.5 0 1.5])
    pause(1)
end
