%% Laboratory notes: Simulation of sampled-data systems in Simulink®
% Riccardo Antonello (riccardo.antonello@unipd.it) and Giulia Michieletto (giulia.michieletto@unipd.it)
% March 10, 2024
% Dept. of Management and Engineering, University of Padova, Italy
%
%
% This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International license.
%
%% Contents
% Introduction
% Continuous-time speed PI control loop
% Sampled-data speed PI control loop
%
%% 1.  Introduction
% Consider the following simplified model of a DC motor:
% (1)   P(s) = \frac{Y(s)}{U_a(s)} = \frac{k}{T\,s+1}\,, \;\quad k = 8.3, \quad T = 0.028
% with the armature voltage u_a\,[\mathrm{V}] as input, and the shaft speed y\,[\mathrm{rad/s}] as output.
% Want to simulate the response of a speed proportional-integral (PI) control loop, assuming to measure the shaft speed (e.g. by using a tachometer).
%
%
%% 2.  Continuous-time speed PI control loop
% Consider first a continuous-time speed PI controller with transfer function:
% (2)   C(s) = \frac{U(s)}{E(s)} = K_P + \frac{K_I}{s}\,,\qquad
% K_P = 8.4\cdot 10^{-3}\,,\quad
% K_I = 1.05
% whose input is the speed error e\,\mathrm{[rpm]}, and the output is the voltage command u\,\mathrm{[V]} provided to the motor voltage driver.
% The motor voltage driver is modelled as a static gain equal to:
% (3)   \frac{U_a(s)}{U(s)} = k_\mathrm{drv} = 0.6
%
% Simulink model of the continuous-time speed PI control loop.

%   conversion gains
rads2rpm = 60/2/pi;
rpm2rads = 2*pi/60;

%   DC motor data
mot.k = 8.3;                %   static gain
mot.T = 0.028;              %   dominant time constant

%   voltage driver data
drv.kdc = 0.6;              %   attenuation gain

%   speed PI controller data
wctrl.kp = 8.4e-3;          %   proportional gain
wctrl.ki = 1.05;            %   integral gain

%   open simulink model   
open_system('speed_ctrl_1.slx');

%
% Note: it is a good practice to organize the model parameters into data structures (e.g. mot, wctrl, etc.), in order to maintain a neat workspace, and avoid conflicts with variable names (e.g. name duplicates).
%
% Simulate the control system response to a step reference input of 300\,\mathrm{rpm}(rotations-per-minute), applied at t=0\,\mathrm{s}.

%   set simulation params
set_param('speed_ctrl_1', ...
    'SolverType', 'Variable-step', ...
    'Solver', 'ode45', ...
    'MaxStep', '0.005', ...
    'StopTime', '0.5');

%   run simulation
sim('speed_ctrl_1');

%
% Plot simualtions results:

%   extract simulation results
t = simres1.time;                           %   time [s]
wm_meas = simres1.signals(1).values(:,1);   %   measured motor speed [rad/s]
wm_ref = simres1.signals(1).values(:,2);    %   motor speed reference [rad/s]
u = simres1.signals(2).values(:,1);         %   driver command [V]

%   plot control response
figure;
subplot(2, 1, 1);
plot(t, wm_meas, 'b-', 'LineWidth', 1.5);
hold on;
plot(t, wm_ref, 'k--', 'LineWidth', 1.5);
grid on;
xlim([0, t(end)]);
ylabel('Speed [rpm]');

subplot(2, 1, 2);
plot(t, u, 'b-', 'LineWidth', 1.5);
ylabel('Voltage [V]')
xlabel('Time [s]');
ylim([0, 10]);
xlim([0, t(end)]);
grid on;

%
% The control system is a LTI system, being an interconnection of LTI blocks, as shown in Fig.1.
%
%
% Figure 1: block diagram of the continuous-time speed PI control loop.
%
% Therefore, its response can be equivalently obtained with the Control Systems Toolbox (CST).
%

%  plant tf
sysP = tf(mot.k, [mot.T, 1])

%  controller tf
sysC = tf([wctrl.kp, wctrl.ki], [1, 0])

%  closed-loop tf (from ref to out)
sysT = feedback( sysC * drv.kdc * sysP * rads2rpm, 1 )

%  closed-loop tf (from ref to controller out)
sysCS = feedback( sysC, drv.kdc * sysP * rads2rpm )

%

%  step response
opt = stepDataOptions('InputOffset', 0, 'StepAmplitude', 300);
t = linspace(0, 0.5, 200);
y = step(sysT, t, opt);
u = step(sysCS, t, opt);

%  plot results
figure;
subplot(2, 1, 1);
plot(t, y, 'b-', 'LineWidth', 1.5);
grid on;
xlim([0, t(end)]);
ylabel('Speed [rpm]');

subplot(2, 1, 2);
plot(t, u, 'b-', 'LineWidth', 1.5);
ylabel('Voltage [V]')
xlabel('Time [s]');
ylim([0, 10]);
xlim([0, t(end)]);
grid on;

%
%% 3.  Sampled-data speed PI control loop
% Consider now the discrete-time speed PI controller with transfer function:
% (4)   C(z) = \frac{U(z)}{E(z)} = K_P + \frac{K_I\,T_s}{z-1}\,,\qquad
% K_P = 8.4\cdot 10^{-3}\,,\quad
% K_I = 1.05\,,\quad
% T_s = 0.01\,\mathrm{s}
% obtained by discretizing the continuous-time controller (2) with the Forward Euler method.

%   sampling time
Ts = 0.01;

%   open simulink model   
open_system('speed_ctrl_2.slx');

%
%
% The resulting system is sampled-data control loop, with a continuous-time plant and a discrete-time speed PI controller (hybrid-time system).
% Note that the Zero-Order Hold block is used either as a discrete-to-continuous interpolator (expected function, as the name suggests ), and as a ideal sampler (counterintuitive). In Simulink, there is no specific blocks to represent the ideal sampler, and one has to resort to the Zero-Order Hold block to perform a continuous-to-discrete transition.
%
% Simulate the control system response to a step reference input of 300\,\mathrm{rpm}, applied at t=0\,\mathrm{s}.

%   set simulation params
set_param('speed_ctrl_2', ...
    'SolverType', 'Variable-step', ...
    'Solver', 'ode45', ...
    'MaxStep', '0.001', ...
    'StopTime', '0.5');

%   run simulation
sim('speed_ctrl_2');

%
% Plot simualtions results:

%   extract simulation results (discrete-time data)
t = simres2.time;                           %   time vector [s]
wm_meas = simres2.signals(1).values(:,1);   %   measured motor speed [rad/s]
wm_ref = simres2.signals(1).values(:,2);    %   motor speed reference [rad/s]
u = simres2.signals(2).values(:,1);         %   driver command [V]

%   extract simulation results (continuous-time data)
tc = simres2c.time;
wm = simres2c.signals.values(:,1);           %   actual motor speed [rad/s]

%   plot control response
figure;
subplot(2, 1, 1);
stairs(t, wm_meas, 'b-', 'LineWidth', 1.5);
hold on;
plot(tc, wm, 'r-', 'LineWidth', 1.5);
stairs(t, wm_ref, 'k--', 'LineWidth', 1.5);
grid on;
xlim([0, 0.5]);
ylabel('Speed [rpm]');
legend('Actual', 'Measured');

subplot(2, 1, 2);
stairs(t, u, 'b-', 'LineWidth', 1.5);
ylabel('Voltage [V]')
xlabel('Time [s]');
ylim([0, 10]);
xlim([0, 0.5]);
grid on;

% Note: being a hybrid-time system, the response of the sampled-data control loop reported above cannot be obtained with the CST.