echo off
% MILLDEMO demonstrates the LQG design tools. The goal is to regulate
% the thickness of a steel beam shaped by a rolling mill stand.

%   Copyright (c) 1986-1999 The Mathworks, Inc. All Rights Reserved.
%   $Revision: 1.6 $  $Date: 1999/01/05 15:21:37 $
clc
clf
randn('seed',24);
help milldemo
echo on
pause % Press any key to continue ...

% The x- and y-axis models of the beam rolling mill includes the 
% following components (see Control System Toolbox User's Guide
% for details):
%
% Hydraulic actuators (with input names "u-x" and "u-y")

Hx = tf(2.4e8,[1 72 90^2],'inputname','u-x')      % x axis 
Hy = tf(7.8e8,[1 71 88^2],'inputname','u-y');     % y axis

% The input thickness/hardness disturbances on each axis are 
% modeled by low-pass filters driven by white noise

Fix = tf(1e4,[1 0.05],'inputname','w-ix');      % x axis 
Fiy = tf(2e4,[1 0.05],'inputname','w-iy');      % y axis

% The rolling cylinder excentricity causes variations in the output
% beam thickness.  These disturbances are modeled by bandpass filters 
% driven by white noise (the filter natural frequency depends on the 
% cylinder revolution period)

Fex = tf([3e4 0],[1 0.125 6^2],'inputn','w-ex');      % x axis
Fey = tf([1e5 0],[1 0.19 9.4^2],'inputn','w-ey');     % y axis

% Gains from force to thickness gap (function of the mill stretch modulus)

gx = 1e-6;
gy = 0.5e-6;


pause % Press any key to continue ...
clc
% LQG design for the x axis:
%
% First build the open-loop model Px.  The plant inputs are the 
% command and the input and excentricity disturbances, and the 
% plant outputs are the deviation "x-gap" to nominal thickness 
% and the measured variation "x-force" of the rolling force 

Px = append([ss(Hx) Fex],Fix);
Px = [-gx gx;1 1] * Px;

% Set plant output names (helps keep track of I/O channels)
set(Px,'outputn',{'x-gap' 'x-force'})


pause % Press any key to continue ...

Px,


pause % Press any key to continue ...
clc
% Perform the state-feedback LQ design.  Penalize low frequency
% content of the thickness gap by adding the low pass filter 
% 30/(s+30) in series with output "x-gap"

lpf = tf(30,[1 30]);
Pxdes = append(lpf,1) * Px;
set(Pxdes,'outputn',{'x-gap*' 'x-force'})

% Design the state-feedback gain using constant weight 1e-4 for 
% command.  Note: keep only the first input "u-x" (command) and 
% first output "x-gap" for design with LQRY

kx = lqry(Pxdes(1,1),1,1e-4)


pause % Press any key to continue ...

% Design the Kalman estimator.  Set measurement noise covariance 
% to 1000 to avoid high-frequency gain.  Keep only the measured 
% output "x-force" (second output) for the estimator design

estx = kalman(Pxdes(2,:),eye(2),1000);


pause % Press any key to continue ...

% Form resulting LQG regulator for x axis and plot its Bode response

Regx = lqgreg(estx,kx);
bode(Regx,{0.1 1e3})

% Note:
% * 0 degree phase in low frequency -> Regulator applies positive 
%   feedback (increased hydraulic force) to reduce extra thickness
%   due to increase in input thickness/hardness (low-freq. disturbance)
% * -180 degree phase near excentricity frequency (6 rad/sec) ->
%   When rolling force increases due to excentricity and beam thickness
%   diminishes as a result, regulator applies negative feedback 
%   (decreased hydraulic force) to offset the loss in beam thickness.


pause % Press any key to continue ...
clc
% Compare the open- and closed-loop responses.
%
% First closed the loop with FEEDBACK (Note: positive feedback required).
% To know which inputs and outputs to connect, look at the I/O names of 
% the plant and regulator

Px.inputn
Regx.outputn

Px.outputn
Regx.inputn

% This indicates that you should connect the first input and second 
% output of Px to the regulator:

clx = feedback(Px,Regx,1,2,+1);

pause % Press any key to continue ...
clc
% Compare the open- and closed-loop Bode responses from input 
% disturbance to output thickness gap (red = open-loop)

bode(Px(1,3),'r',clx(1,3),'b',{0.01 100})

% The LQG regulator reduces the low frequency gain by about 20 dB 


pause % Press any key to continue ...

% Compare the open- and closed-loop Bode responses from excentricity 
% disturbance to output thickness gap (red = open-loop)

bode(Px(1,2),'r',clx(1,2),'b',{1,100})

% The LQG regulator attenuates excentricity resonance by about 20 dB 


pause % Press any key to continue ...
clc
% Now compare open- and closed-loop time responses to unit-covariance
% white noise disturbances w-ex and w-ix.  Sample the noise input every
% dt=0.01 second.  The equivalent discrete noise covariance is then 1/dt

dt = 0.01;
if isstudent, dt = 0.1; end
t = 0:dt:50;                             % input time samples
wx = sqrt(1/dt) * randn(2,length(t));    % sampled driving noise

% Simulate response with LSIM (green = open-loop)

lsim(Px(1,2:3),'g',clx(1,2:3),'b',wx,t)


pause % Press any key to continue ...
clc
% Perform a similar design for the y axis. Form the open-loop model:

Py = append([ss(Hy) Fey],Fiy);
Py = [-gy gy;1 1] * Py;
set(Py,'outputn',{'y-gap' 'y-force'})

% State-feedback gain design:

Pydes = append(lpf,1) * Py;
set(Pydes,'outputn',{'y-gap*' 'y-force'})
ky = lqry(Pydes(1,1),1,1e-4);

% Kalman estimator design:

esty = kalman(Pydes(2,:),eye(2),1000);

% Form LQG regulator and closed-loop:
Regy = lqgreg(esty,ky);
cly = feedback(Py,Regy,1,2,+1);

pause % Press any key to continue ...
clc
% Compare the open- and closed-loop time responses to white noise 
% disturbances w-ey and w-iy using LSIM (green = open-loop)

dt = 0.01;
if isstudent, dt = 0.1; end
t = 0:dt:50;
wy = sqrt(1/dt) * randn(2,length(t));

lsim(Py(1,2:3),'g',cly(1,2:3),'b',wy,t)


pause % Press any key to continue ...
clc
% You have designed SISO regulators for each axis independently.
% In the two-dimensional rolling mill, however, there is some
% cross-coupling between the rolling forces along each axis (see
% Control System Toolbox User's Guide for details).  The coupling
% coefficients are:

gxy = 0.1; 
gyx = 0.4; 


pause % Press any key to continue ...

% Build the two-axis open-loop model P with cross-coupling terms:

P = append(Px,Py);                % Append x- and y-axis models
P = P([1 3 2 4],[1 4 2 3 5 6]);   % Reorder inputs and outputs 

P.inputn
P.outputn

% The coupling terms amount to a linear combination of the outputs 
% of Px and Py (x- and y-axis models).  Define and add the 
% corresponding static gain

CouplingGain = [eye(2) [0 gyx*gx;gxy*gy 0] ; zeros(2) [1 -gyx;-gxy 1]];
Pc = CouplingGain * P;
Pc.outputn = P.outputn;


pause % Press any key to continue ...
clc
% Using this detailed two-axis model, perform open- and closed-loop 
% simulations of the x/y thickness gaps.  Use the two SISO regulators 
% in closed-loop simulation.

% Close the two SISO feedback loops

cl = feedback(Pc,append(Regx,Regy),1:2,3:4,+1);

% Compare open- and closed-loop responses of x/y axes 
% (green = open loop)

wxy = [wx ; wy];
lsim(Pc(1:2,3:6),'g',cl(1:2,3:6),'b',wxy,t)

% The coupling between axes severely deteriorates the regulation 
% performance along the x axis.  Designing for one loop at a time
% is therefore inadequate for this level of cross-coupling.



pause % Press any key to continue ...
clc
% To overcome the cross-coupling effects, design a MIMO LQG regulator 
% for the complete x/y-axis model P.

% Add low-pass filter used to penalize x/y-gap in LQRY design

Pdes = append(lpf,lpf,eye(2)) * Pc;
Pdes.outputn = Pc.outputn;

% Design the LQ gain and Kalman estimator using
%   * x-force and y-force as measurements
%   * u-x and u-y as commands
%   * Q = I , R = 1e-4*I

k = lqry(Pdes(1:2,1:2),eye(2),1e-4*eye(2))

est = kalman(Pdes(3:4,:),eye(4),1000*eye(2));

RegMIMO = lqgreg(est,k);



pause % Press any key to continue ...
clc
% Plot the singular value response of this two-input/two-output
% regulator

sigma(RegMIMO)

pause % Press any key to continue ...


% Joint-axis simulation with MIMO regulator and cross-coupling effects.
% Compare open- and closed-loop time responses (green = open loop):

cl = feedback(Pc,RegMIMO,1:2,3:4,+1);

lsim(Pc(1:2,3:6),'g',cl(1:2,3:6),'b',wxy,t)

% The MIMO design is a clear improvement over the separate SISO designs 
% for each axis. In particular, the level of x/y thickness variation 
% is now comparable to that obtained in the decoupled case. 
% This example illustrates the benefits of direct MIMO design for 
% multivariable systems.

echo off
