Zakir Husain Delhi College

University of Delhi

Control Systems

Submitted to: Dr. Manish Kumar

Submitted By: Tanmoi Saikia

Roll: 22/1200

Submitted By: Tanmoi Saikia

 INDEX

Sl. Experiment Date Sign

No
1 To study response of systems for
various standard test input signals.
2 To study position and speed control
of DC motor.
3 To study time and frequency domain
specifications of a control system.
4 To plot Bode, Root locus and
Nyquist plots and determine
stability.
5 To study the effect of PI, PD and
PID controller on closed loop
systems
6 State space analysis for a given
Transfer function
 Experiment 1

Aim: To study response of systems for various standard test input signals.

(a) Step Response of First-Order Systems for Different Settling Times

Source Code:

clc;
clear;
close;

// Time vectors
t1 = 0:0.01:10; // For Ts = 4s
t2 = 0:0.001:1; // For Ts = 0.4s
t3 = 0:0.1:100; // For Ts = 40s

// Given R values
R1 = 2; R2 = 1; R3 = 4;

// Desired settling times

Ts1 = 4; Ts2 = 0.4; Ts3 = 40;

// Compute tau and C for each

tau1 = Ts1 / 4; tau2 = Ts2 / 4; tau3 = Ts3 / 4;
C1 = tau1 / R1; C2 = tau2 / R2; C3 = tau3 / R3;

// Compute alpha
alpha1 = 1 / tau1; alpha2 = 1 / tau2; alpha3 = 1 / tau3;

// Define Laplace variable

s = poly(0,'s');

// Transfer functions
G1 = syslin('c', alpha1/(s + alpha1));
G2 = syslin('c', alpha2/(s + alpha2));
G3 = syslin('c', alpha3/(s + alpha3));

// Step responses
y1 = csim('step', t1, G1);
y2 = csim('step', t2, G2);
y3 = csim('step', t3, G3);

// Plot
clf;
subplot(3,1,1);
plot(t1, y1, 'r');
title("Step Response (R="+string(R1)+"Ω, C="+string(C1)+"F)");
xlabel('Time (s)'); ylabel('Output'); xgrid();

Submitted By: Tanmoi Saikia

subplot(3,1,2);
plot(t2, y2, 'g');
title("Step Response (R="+string(R2)+"Ω, C="+string(C2)+"F)");
xlabel('Time (s)'); ylabel('Output'); xgrid();

subplot(3,1,3);
plot(t3, y3, 'b');
title("Step Response (R="+string(R3)+"Ω, C="+string(C3)+"F)");
xlabel('Time (s)'); ylabel('Output'); xgrid();

Output:

(b) Time Response Analysis of a First-Order System for Standard Test

Signals

Source Code:
// Clear environment
clc;
clear;
close;

// Time vector
t = 0:0.01:10; // Time from 0 to 10 seconds
dt = t(2) - t(1); // Time step = 0.01

// Define Transfer Function: H(s) = 1 / (s + 1)

s = poly(0, 's');
num = 1;
den = s + 1;
sys = syslin('c', num, den); // Continuous-time system

Submitted By: Tanmoi Saikia

// 1. Step Input
u_step = ones(t); // Step input
y_step = csim(u_step, t, sys); // System response
e_step = u_step - y_step; // Error function

// 2. Ramp Input
u_ramp = t; // Ramp input
y_ramp = csim(u_ramp, t, sys); // System response
e_ramp = u_ramp - y_ramp; // Error function

// 3. Impulse Input (approximate)

u_impulse = zeros(t);
u_impulse(1) = 1 / dt; // Approximate unit impulse
y_impulse = csim(u_impulse, t, sys); // System response
e_impulse = u_impulse - y_impulse; // Error function

// Plotting Inputs, Outputs, and Errors

scf(0); // Create new figure

//-------- STEP --------

subplot(3,3,1); // Step Input
plot(t, u_step, 'b');
title('Step Input');
xlabel('Time (s)');
ylabel('Amplitude');
xgrid();

subplot(3,3,2); // Step Output

plot(t, y_step, 'r');
title('System Response to Step');
xlabel('Time (s)');
ylabel('Output');
xgrid();

subplot(3,3,3); // Step Error

plot(t, e_step, 'm');
title('Error e(t) for Step');
xlabel('Time (s)');
ylabel('Error');
xgrid();

//-------- RAMP --------

subplot(3,3,4); // Ramp Input
plot(t, u_ramp, 'b');
title('Ramp Input');
xlabel('Time (s)');
ylabel('Amplitude');
xgrid();

subplot(3,3,5); // Ramp Output

plot(t, y_ramp, 'r');
title('System Response to Ramp');
xlabel('Time (s)');
ylabel('Output');

Submitted By: Tanmoi Saikia

xgrid();

subplot(3,3,6); // Ramp Error

plot(t, e_ramp, 'm');
title('Error e(t) for Ramp');
xlabel('Time (s)');
ylabel('Error');
xgrid();

//-------- IMPULSE --------

subplot(3,3,7); // Impulse Input
plot(t, u_impulse, 'b');
title('Impulse Input (approx)');
xlabel('Time (s)');
ylabel('Amplitude');
xgrid();

subplot(3,3,8); // Impulse Output

plot(t, y_impulse, 'r');
title('System Response to Impulse');
xlabel('Time (s)');
ylabel('Output');
xgrid();

subplot(3,3,9); // Impulse Error

plot(t, e_impulse, 'm');
title('Error e(t) for Impulse');
xlabel('Time (s)');
ylabel('Error');
xgrid();

Output:
For Step Signal

Submitted By: Tanmoi Saikia

For Ramp Signal

For Impulse Signal

Submitted By: Tanmoi Saikia

For 2nd Order Systems
(i) Unit Step Input
Source Code:
clear; clf;
wn = 2*%pi; t = linspace(0,3,5000);
ub = 1.02; lb = 0.98;

function plot_response(zeta, pos)

subplot(2,2,pos);
select zeta
case 0 then
c = 1 - cos(wn*t); ttl = "Undamped (ζ=0)";
case 1 then
c = 1 - (1 + wn*t).*exp(-wn*t); ttl = "Critically Damped (ζ=1)";
else
if zeta<1 then
wd = wn*sqrt(1-zeta^2); th = atan(sqrt(1-zeta^2)/zeta);
c = 1 - exp(-zeta*wn*t).*sin(wd*t+th)/sqrt(1-zeta^2);
ttl = "Underdamped (ζ="+string(zeta)+")";
else
s1 = -zeta*wn + wn*sqrt(zeta^2-1);
s2 = -zeta*wn - wn*sqrt(zeta^2-1);
c = 1 - (s2*exp(s1*t)-s1*exp(s2*t))/(s2-s1);
ttl = "Overdamped (ζ="+string(zeta)+")";
end
end

plot(t,c,'k',t,ones(t),'k--');
title(ttl); xlabel("Time (s)"); ylabel("Amplitude"); xgrid(1);

idx90 = find(c>=0.9,1); tr = %inf; if ~isempty(idx90) then tr=t(idx90); plot(tr,c(idx90),'ko'); end

idxs = find((c>=lb)&(c<=ub)); ts = %inf; if ~isempty(idxs) then ts=t(idxs(1)); plot([ts ts],[lb ub],'k'); end
[cmax,imax]=max(c); Mp=(cmax-1)*100; plot(t(imax),cmax,'ko');
 if zeta==0 then info=["Period="+msprintf("%.2f",2*%pi/wn);"M_p=100%";"No Settling"];
elseif zeta<1 then
info=["t_r="+msprintf("%.2f",tr);"t_p="+msprintf("%.2f",t(imax));"t_s="+msprintf("%.2f",ts);"M_p="+msprintf("%.
1f%%",Mp)];
else info=["t_r="+msprintf("%.2f",tr);"t_s="+msprintf("%.2f",ts);"M_p=0%"]; end
xstring(1.5,0.6,info);
endfunction

plot_response(0.3,1);
plot_response(1,2);
plot_response(1.5,3);
plot_response(0,4);

Output:

(ii) Unit Impulse Input

Source Code:
clear; clf;
// Common parameters
wn = 2 * %pi; // Natural frequency (rad/s)
t = linspace(0, 3, 5000); // Time vector
// --- Function to compute & plot impulse responses ---
function plot_impulse(zeta, pos)
subplot(2, 2, pos);

select zeta
case 0 then
c = wn * sin(wn * t);
ttl = "Undamped (ζ=0)";
case 1 then
c = wn^2 * t .* exp(-wn * t);
ttl = "Critically Damped (ζ=1)";
else
if zeta < 1 then
wd = wn * sqrt(1 - zeta^2);
c = (wn / sqrt(1 - zeta^2)) * exp(-zeta * wn * t) .* sin(wd * t);
ttl = "Underdamped (ζ=" + string(zeta) + ")";
else
s1 = -zeta*wn + wn*sqrt(zeta^2 - 1);
s2 = -zeta*wn - wn*sqrt(zeta^2 - 1);
c = (wn^2 / (s2 - s1)) * (exp(s1*t) - exp(s2*t));
ttl = "Overdamped (ζ=" + string(zeta) + ")";
end
end

// Plot response
h = plot(t, c);
[Link] = 1; [Link] = 2.5;
xgrid(1);
xlabel("Time (s)");
ylabel("Amplitude");
title(ttl);
endfunction

// --- Run for all damping cases ---

plot_impulse(0.3,1); // Underdamped
plot_impulse(1,2); // Critically damped
plot_impulse(1.5,3); // Overdamped
plot_impulse(0,4); // Undamped

Output:

(iii) Ramp Input

Source Code:
clear; clf;

// Common parameters
wn = 2 * %pi; // Natural frequency (rad/s)
t = linspace(0, 5, 5000); // Time vector
y_ramp = t; // Ideal ramp (input reference)

// --- Function to compute & plot ramp response ---

function plot_ramp(zeta, pos)
subplot(2, 2, pos);

select zeta
case 0 then
// Undamped: infinite oscillations about ramp line
c = t - sin(wn * t) / wn;
ttl = "Undamped (ζ=0)";
case 1 then
// Critically damped
c = t - exp(-wn * t) .* (t + 1/wn);
ttl = "Critically Damped (ζ=1)";
else
if zeta < 1 then
wd = wn * sqrt(1 - zeta^2);
phi = atan(sqrt(1 - zeta^2) / zeta);
// Ramp response for underdamped case
c = t - (exp(-zeta * wn * t) / (wn * sqrt(1 - zeta^2))) .* ...
(sin(wd * t + phi));
ttl = "Underdamped (ζ=" + string(zeta) + ")";
else
s1 = -zeta*wn + wn*sqrt(zeta^2 - 1);
s2 = -zeta*wn - wn*sqrt(zeta^2 - 1);
c = t - (exp(s1*t)/s1 - exp(s2*t)/s2)/(s2 - s1)*wn^2;
ttl = "Overdamped (ζ=" + string(zeta) + ")";
end
end

// Plot actual and reference ramp

h = plot(t, c, "k"); [Link] = 2.5;
plot(t, y_ramp, 'r--'); // Reference ramp in red dashed
xgrid(1);
xlabel("Time (s)");
ylabel("Amplitude");
title(ttl);
 // Steady-state error
ess = 2 * zeta / wn;
xstring(2, 0.7 * max(c), "e_ss = " + msprintf("%.3f", ess));
endfunction

// --- Run all damping cases ---

plot_ramp(0.3,1); // Underdamped
plot_ramp(1,2); // Critically damped
plot_ramp(1.5,3); // Overdamped
plot_ramp(0,4); // Undamped

Output:
(iv) Sinusoidal Input
Source Code:
clear; clf;

// Common parameters
wn = 2 * %pi; // Natural frequency (rad/s)
t = linspace(0, 5, 5000); // Time vector (0–5 s)
w_input = 1.5 * %pi; // Input signal frequency (rad/s)
r = sin(w_input * t); // Input sinusoid

// --- Function to compute & plot sinusoidal response ---

function plot_sine(zeta, pos)
subplot(2, 2, pos);

// Magnitude and phase of frequency response

mag = wn^2 / sqrt((wn^2 - w_input^2)^2 + (2*zeta*wn*w_input)^2);
phi = -atan((2*zeta*wn*w_input) / (wn^2 - w_input^2));

// Steady-state output
c = mag * sin(w_input * t + phi);

// Plot input and output

plot(t, r, 'r--'); // input in red dashed
h = plot(t, c, 'k'); [Link] = 2.5; // output in black

xgrid(1);
xlabel("Time (s)");
ylabel("Amplitude");
ttl = "Sinusoidal Input Response (ζ=" + string(zeta) + ")";
title(ttl);

// Gain and phase info

gain_db = 20 * log10(mag);
info = ["|T(jω)| = " + msprintf("%.3f", mag);
"φ = " + msprintf("%.2f°", phi * 180 / %pi);
 "Gain = " + msprintf("%.2f dB", gain_db)];
xstring(1, 0.7, info);
endfunction

// --- Run cases ---

plot_sine(0,1); // Undamped
plot_sine(0.3,2); // Underdamped
plot_sine(1,3); // Critically damped
plot_sine(1.5,4); // Overdamped

Output:
 Experiment 2

Aim: To study the characteristics of a DC Motor

Source Code:
clc;
clear;
close;
Vm = 0:0.5:12; // Motor voltage (V)
Km = 348.25; // Motor constant (rpm/volt)
N = Km * Vm; // Linear region
N(N > 2300) = 2300; // Saturation near 2300 rpm

scf(1);
plot(Vm, N, 'b-o', 'LineWidth', 2);
xgrid(1);
xlabel("Motor Voltage V_m (V)", "fontsize", 3);
ylabel("Speed N (rpm)", "fontsize", 3);
title("Motor Characteristics", "fontsize", 4);
set(gca(), "data_bounds", [0, 0; 12, 3000]);

// Annotation
xstring(5, 1500, "K_m = 348.25 rpm/volt");
xstring(5, 1300, "= 36.458 rad/volt·s");

N2 = 0:500:3000; // Speed (rpm)

Kt = 1/1825; // Tacho constant (volt/RPM)
Vt = Kt * N2; // Linear tacho voltage

scf(2);
plot(N2, Vt, 'r-s', 'LineWidth', 2);
xgrid(1);
xlabel("Speed N (rpm)", "fontsize", 3);
ylabel("Tacho Voltage V_t (V RMS)", "fontsize", 3);
title("Tacho Characteristics", "fontsize", 4);
set(gca(), "data_bounds", [0, 0; 3000, 1.6]);

// Annotation
xstring(1000, 1, "K_t = 1/1825 volt/RPM");
xstring(1000, 0.9, "= 0.00523 volt·s/rad");

Output:
 Experiment 3

Aim: To study time and frequency domain specifications of a control

system.
Source Code:
clc;
clear;
close;

// -------------------- Parameters --------------------

T = 10; // Fundamental period (s)
w0 = 2*%pi / T; // Fundamental angular frequency (rad/s)
N = 40; // Number of harmonics
dt = 1e-3; // integration step

// Time vector
t = -5:0.01:5;
Nt = length(t);

// Integration vector for one period

t_integ = 0:dt:(T-dt);
f = double((t_integ >= 0) & (t_integ < T/2)); // one-period unit step

// -------------------- Fourier Coefficients --------------------

a0 = (2/T) * sum(f) * dt;
an = zeros(1, N);
bn = zeros(1, N);

for n = 1:N
cosvec = cos(n*w0*t_integ);
sinvec = sin(n*w0*t_integ);
an(n) = (2/T) * sum(f .* cosvec) * dt;
bn(n) = (2/T) * sum(f .* sinvec) * dt;
end
// -------------------- Reconstruction --------------------
x_reconstructed = (a0/2) * ones(1, Nt);
for n = 1:N
x_reconstructed = x_reconstructed + an(n) * cos(n*w0*t) + bn(n) * sin(n*w0*t);
end

u = double(t >= 0);

// -------------------- Time-domain Plot --------------------

clf();
subplot(2,2,1);
plot(t, u, 'r', 'LineWidth', 3);
title('Original Unit Step (u(t))');
xlabel('Time (s)');
ylabel('Amplitude');
xgrid();

subplot(2,2,2);
plot(t, x_reconstructed, 'b', 'LineWidth', 3);
title('Reconstructed using Fourier Series');
xlabel('Time (s)');
ylabel('Amplitude');
xgrid();

// -------------------- Frequency-domain (Magnitude + Phase) --------------------

n_vec = 1:N;
freq = n_vec * w0 / (2*%pi);
magnitude = sqrt(an.^2 + bn.^2);
phase = atan(-bn ./ an); // phase of Cn = an - j*bn

// Magnitude Spectrum
subplot(2,2,3);
plot(freq, magnitude, 'o-', 'LineWidth', 2);
title('Magnitude Spectrum |C_n|');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
xgrid();

// Phase Spectrum
subplot(2,2,4);
plot(freq, phase, 's-', 'LineWidth', 2);
title('Phase Spectrum ∠C_n');
xlabel('Frequency (Hz)');
ylabel('Phase (rad)');
xgrid();

// -------------------- Info --------------------

disp("------ Fourier series summary ------");
disp("Fundamental frequency (Hz) : " + string(w0/(2*%pi)));
disp("a0/2 (DC value in reconstruction): " + string(a0/2));

Output:
 Experiment 4

Aim: To plot Bode, Root locus and Nyquist plots and determine stability.
(I) For Bode Plot
s = %s;
G = 1/(s+1); // Transfer function
w = logspace(-2,2,200); // Frequency range
jw = %i*w; // s = jw
Gjw = horner(G, jw); // Evaluate G(jw)
magnitude_dB = 20*log10(abs(Gjw));
phase_deg = atand(imag(Gjw), real(Gjw));

subplot(2,1,1);
plot2d("ln", w, magnitude_dB);
xgrid(); title('Bode Magnitude Plot');
ylabel('Magnitude (dB)');
subplot(2,1,2);
plot2d("ln", w, phase_deg);
xgrid(); title('Bode Phase Plot');
ylabel('Phase (deg)'); xlabel('Frequency (rad/s)');

Output:
 (ii) For Nyquist Plot
s = %s;
G = 1/(s+1); // Transfer function
w = logspace(-2,2,200); // Frequency range
jw = %i*w; // s = jw
Gjw = horner(G, jw); // Evaluate G(jw)
// Nyquist plot (Real vs Imaginary)
plot(real(Gjw), imag(Gjw), 'b-', 'LineWidth', 2);
plot(real(Gjw($:-1:1)), -imag(Gjw($:-1:1)), 'b--', 'LineWidth', 1); // Negative frequencies
xgrid(); title('Nyquist Plot');
xlabel('Real'); ylabel('Imaginary');
axis equal;

Output:
 (iii) For Root Locus Plot

s = %s;
G = 1/(s*(s+1)); // Transfer function
k = logspace(-2,2,100); // Gain values

clf();
for i = 1:length(k)
char_eq = %s^2 + %s + k(i); // Characteristic equation: 1 + k*G(s) = 0
roots_val = roots(char_eq); // Find roots
plot(real(roots_val), imag(roots_val), 'b.');
end
xgrid(); title('Root Locus Plot');
xlabel('Real Axis'); ylabel('Imaginary Axis');
axis equal;

Output:
 Experiment 5

Aim: To study the effect of PI, PD and PID controller on closed loop
systems
Source Code:
clc; clear; close;

// Laplace variable and plant

s = poly(0, 's');
G = syslin('c', 1/s);
t = 0:0.01:10;
K_values = [0.5, 1, 2, 5];

// Controller parameters
Ki = 1; Kd = 0.5;

// Controller types and equations

controllers = ["PI", "PD", "PID"];
C_expr = ["Kp + Ki/s", "Kp + Kd*s", "Kp + Ki/s + Kd*s"];

for idx = 1:3

scf(idx); clf(); // Open a new figure window
legend_entries = [];

for i = 1:length(K_values)
Kp = K_values(i);
execstr("C = " + C_expr(idx) + ";");
T = syslin('c', (C*G)/(1 + C*G));
y = csim('step', t, T);
plot(t, y, 'LineWidth', 2);
legend_entries($+1) = "K=" + string(Kp);
if i == 1 then
title(controllers(idx) + " Controller Responses for Different K");
xlabel("Time (s)");
 ylabel("Output");
xgrid();
end
end
legend(legend_entries, "in_lower_right");
end

Output:
 Experiment 6

Aim: State space analysis for a given Transfer function

Source Code:
clc;
clear;
close;
t = -2:0.01:2; // includes negative and positive values
x1 = -1 + 2*exp(t); // x₁(t)
x2 = 2 .* t .* exp(t); // x₂(t)
y = -1 + 2*exp(t); // y(t)

scf(1);
plot(t, x1, 'r', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('x₁(t)');
title('State Variable: x₁(t)');
xgrid();
legend('x₁(t)', 'location', 'topleft');

// Mark point at t = 0
x1_0 = -1 + 2*exp(0);
plot(0, x1_0, 'ko');
xstring(0, x1_0, " (0,1)");
xstring(-1.8, 5, "x₁(t) = -1 + 2e^t");

scf(2);
plot(t, x2, 'b', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('x₂(t)');
title('State Variable: x₂(t)');
xgrid();
legend('x₂(t)', 'location', 'topleft');
xstring(-1.8, 10, "x₂(t) = 2te^t");
scf(3);
plot(t, y, 'k', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('y(t)');
title('Output: y(t)');
xgrid();
legend('y(t)', 'location', 'topleft');

y_0 = -1 + 2*exp(0);
plot(0, y_0, 'ko');
xstring(0, y_0, " (0,1)");

xstring(-1.8, 5, "y(t) = -1 + 2e^t");

mprintf("\nAnalytical Solutions:\n");
mprintf("x₁(t) = -1 + 2e^t\n");
mprintf("x₂(t) = 2te^t\n");
mprintf("y(t) = -1 + 2e^t\n");
mprintf("\nAt t = 0: x₁(0) = %.2f, x₂(0) = %.2f, y(0) = %.2f\n", x1_0, x2(find(t==0)), y_0);

Output: