Course Name:

Control System
Submitted by:
FARRUKH MEHMOOD (2022-EE-018)

Submitted to:
Dr. Habib Wajid
Section :
A
Lab # 7
Task 1:
Code:
%% Task 1: Magnetically Suspended Ball - System Definition
A = [ 0 1 0;
980 0 -2.8;
0 0 -100];
B = [ 0;
0;
100];
C = [1 0 0];
D = 0;
sys_open_loop = ss(A, B, C, D);
t = 0:0.01:2;
x0 = [0.01; 0; 0];
u = zeros(size(t));
% Task 1.a: Open-loop System Stability
eigenvalues_A = eig(A);
% Task 1.b: Unstable System Response to Non-Zero Initial Condition
[y_open_loop, t_open_loop, x_open_loop] = lsim(sys_open_loop, u, t, x0);
% Task 1.c: Controllability and Observability
n = size(A, 1); % Number of states
Controllability_Matrix = ctrb(A, B);
rank_Controllability = rank(Controllability_Matrix);
Observability_Matrix = obsv(A, C);
rank_Observability = rank(Observability_Matrix);
% Task 1: Pole Placement Controller Design - First Attempt
p1_first=-10+10i;
p2_first=-10-10i;
p3_first=-50;
K_first_attempt = place(A,B,[p1_first p2_first p3_first]);
sys_cl_first_attempt = ss(A-B*K_first_attempt,B,C,0); % D is 0 as per manual's sys_cl definition
% Using the general t and x0 for simulation
[y_cl_first, t_cl_first, x_cl_first] = lsim(sys_cl_first_attempt, u, t, x0);
% Second set of desired closed-loop poles: -20 +/- 20i, -100.
p1_second=-20+20i;
p2_second=-20-20i;
p3_second=-100;
K_second_attempt = place(A,B,[p1_second p2_second p3_second]);
sys_cl_second_attempt = ss(A-B*K_second_attempt,B,C,0);
% Using the general t and x0 for simulation
[y_cl_second, t_cl_second, x_cl_second] = lsim(sys_cl_second_attempt, u, t, x0);
% Task 1.d: Step Response and Reference Scaling (Nbar)
u_step_val = 0.001;
t_step = 0:0.001:2; % Time vector for step response
u_step = u_step_val * ones(size(t_step)); % Step input (constant input for step response)
x0_step = [0; 0; 0];
[y_step_no_Nbar, t_step_no_Nbar, x_step_no_Nbar] = lsim(sys_cl_second_attempt, u_step, t_step, x0_step);

% Simulate step response with Nbar

[y_step_with_Nbar, t_step_with_Nbar, x_step_with_Nbar] = lsim(sys_cl_second_attempt, Nbar * u_step, t_step,
x0_step);
% Task 1.e: Observer Design and Simulation
op1=-100;
op2=-101;
op3=-102;
P_observer = [op1 op2 op3];
L = place(A', C', P_observer)';
disp(L);
At=[A-B*K_second_attempt B*K_second_attempt; zeros(size(A)) A-L*C];
Bt=[B*Nbar; zeros(size(B))];
Ct=[C zeros(size(C))];
Dt_combined = 0;
sys_obs_combined = ss(At, Bt, Ct, Dt_combined);

t_observer_plot = 0:1E-6:0.1;
x0_observer_plot_manual=[0.01 0.5 -5];
x0_obs_combined_manual = [x0_observer_plot_manual x0_observer_plot_manual];
[y_obs_manual, ~, x_states_error_manual] = lsim(sys_obs_combined, zeros(size(t_observer_plot)),
t_observer_plot, x0_obs_combined_manual);
n_states = 3; % Number of states for the original system
e_manual=x_states_error_manual(:,n_states+1:end); % Estimation error
x_actual_manual=x_states_error_manual(:,1:n_states); % Actual states
x_estimated_manual=x_actual_manual-e_manual; % Estimated states
h=x_actual_manual(:,1); h_dot=x_actual_manual(:,2); i=x_actual_manual(:,3);
h_est=x_estimated_manual(:,1); h_dot_est=x_estimated_manual(:,2); i_est=x_estimated_manual(:,3);
figure;
plot(t_observer_plot, h,'-r',t_observer_plot, h_est,':r',...
t_observer_plot, h_dot,'-b',t_observer_plot,h_dot_est,':b',...
t_observer_plot, i,'-g',t_observer_plot,i_est,':g');
legend ('h', 'h_{est}', 'hdot', 'hdot_{est}', 'i','i_{est}'); % Corrected legend for proper display
xlabel('Time (sec)');
ylabel('State Values');
title('State Variables and Observer Estimates');
grid on;
disp('Plots showing actual and estimated state variables are generated.');
disp('The observer estimates should converge to actual state variables.');

Output:
Task 2:
Code:

% Task 2: Cruise Control - System Definition

m = 1000; % kg (vehicle mass)
b = 50; % N.s/m (damping coefficient)
A = -b/m; % A = -50/1000 = -0.05
B = 1/m; % B = 1/1000 = 0.001
C = 1; % Output is velocity
D = 0; % No direct path
sys_open_loop_cc = ss(A, B, C, D);
t = 0:0.1:20;
x0_step_sim = 0;
% Task 2.a: Pole Placement Controller Design
desired_pole = -1.5;
K_cc = place(A, B, desired_pole); % state feedback gain
A_cl = A - B*K_cc;
sys_cl_cc = ss(A_cl, B, C, D);
u_unit_step_ref = 1 * ones(size(t));
[y_cl_no_Nbar, t_cl_no_Nbar] = lsim(sys_cl_cc, u_unit_step_ref, t, x0_step_sim);
% Task 2.b: Compute Nbar to Eliminate Steady-State Error
Nbar_calculated = 1 / (C * inv(A - B*K_cc) * B);
% Define system with scaled input
sys_cl_with_Nbar = ss(A - B*K_cc, B*Nbar_calculated, C, D);
% Simulate response for 10 m/s desired speed
desired_steady_state_speed = 10;
u_input_with_Nbar = ones(size(t)); % unit step input
[y_cl_with_Nbar, t_cl_with_Nbar] = lsim(sys_cl_with_Nbar, u_input_with_Nbar, t, x0_step_sim);
y_cl_with_Nbar = y_cl_with_Nbar * desired_steady_state_speed;
% Step Response Characteristics
step_characteristics = stepinfo(y_cl_with_Nbar, t);
rise_time = step_characteristics.RiseTime;
overshoot_percent = step_characteristics.Overshoot;
% Steady-state value = last value of y
final_value_10mps = y_cl_with_Nbar(end);
steady_state_error_abs = abs(desired_steady_state_speed - final_value_10mps);
steady_state_error_percent = (steady_state_error_abs / desired_steady_state_speed) * 1;
% Display Results
disp(' ');
disp('--- Design Completion Verification ---');
disp(['Desired Steady-State Speed: ', num2str(desired_steady_state_speed), ' m/s']);
disp(['Actual Steady-State Speed (from simulation): ', num2str(final_value_10mps), ' m/s']);
disp(' ');
disp(['Rise Time: ', num2str(rise_time), ' s (Requirement: < 5s)']);
disp(['Overshoot: ', num2str(overshoot_percent), '% (Requirement: < 10%)']);
disp(['Steady-State Error: ', num2str(steady_state_error_percent), '% (Requirement: < 2%)']);

Output:

Task 3:
Code:
clear; clc;

% Task 3: DC Motor Speed - System Definition

% Given parameters
J = 0.01; % kg·m^2
b = 0.1; % N·m·s
K = 0.01; % motor constant
R = 1; % Ohm
L = 0.5; % Henry
% State-space matrices
A = [-b/J K/J;
-K/L -R/L];
B = [0;
1/L];
C = [1 0]; % Output is angular velocity
D = 0;
sys_motor = ss(A, B, C, D);
% Step 1: Check System Order and Controllability
n_states = size(A, 1)
rank_ctrb = rank(ctrb(A, B))
% Step 2: Desired Pole Placement
% Desired poles: -5 ± j (settling time ≈ 0.8 s, OS < 5%)
p1 = -5 + 1j;
p2 = -5 - 1j;
K_motor = place(A, B, [p1, p2]);
A_cl = A - B*K_motor;
sys_cl = ss(A_cl, B, C, D);
% Step 3: Simulate Step Response (No Nbar)
t = 0:0.01:3;
u = ones(size(t)); % 1 rad/sec reference
x0 = [0; 0];
[y, t, ~] = lsim(sys_cl, u, t, x0);
% Step 4: Compute Nbar to Eliminate Steady-State Error
Nbar = 1 / (C * inv(A - B*K_motor) * B);
sys_cl_scaled = ss(A - B*K_motor, B*Nbar, C, D);
% Simulate response with Nbar
[y_scaled, ~, ~] = lsim(sys_cl_scaled, u, t, x0);
% Step 5: Step Response Characteristics
step_info = stepinfo(y_scaled, t);
final_value = y_scaled(end);
ss_error = abs(0 - final_value) ;
disp('--- Performance Check ---');
disp(['Rise Time: ', num2str(step_info.RiseTime), ' s (Req: <2s)']);
disp(['Overshoot: ', num2str(step_info.Overshoot), '% (Req: <5%)']);
disp(['Steady-State Error: ', num2str(ss_error), '% (Req: <1%)']);

Output:

Task 4:
Code:
% Task 4: DC Motor Position - System Parameters
J = 3.2284e-6;
b = 3.5077e-6;
K = 0.0274;
R = 4;
L = 2.75e-6;
% Step 1: State-Space Representation
A = [0 1 0;
0 -b/J K/J;
0 -K/L -R/L];
B = [0;
0;
1/L];
C = [1 0 0]; % Output is position
D = 0;
sys_motor_pos = ss(A, B, C, D);
% Step 2: Check Order & Controllability
n_states = size(A,1)
rank_ctrb = rank(ctrb(A,B))
% Step 3: Pole Placement (Initial Controller, No Integrator Yet)
p1 = -100 + 100i;
p2 = -100 - 100i;
p3 = -200;
Kc = place(A, B, [p1, p2, p3]);
% Step 4: Augment System with Integrator for Zero Steady-State Error
Aa = [0 C;
zeros(3,1) A];
Ba = [0;
B];
Ca = [0 C];
% Desired poles: add one fast integrator pole
P_aug = [-100+100i, -100-100i, -200, -300]; % 4 poles
Ka = place(Aa, Ba, P_aug); % Combined gain
Ki = Ka(1); % Integrator gain
K = Ka(2:end); % State-feedback gain
% Step 5: Closed-Loop System with Integral Action
A_cl = [A - B*K, B*Ki;
-C, 0];
B_cl = [zeros(3,1);
1];
C_cl = [C, 0];
D_cl = 0;
sys_aug = ss(A_cl, B_cl, C_cl, D_cl);
% Step 6: Performance Evaluation
info = stepinfo(y, t);
final_value = y(end);
ss_error = abs(1 - final_value) * 100;
disp('--- Design Verification ---');
disp(['Settling Time: ', num2str(info.SettlingTime), ' s (Req: < 0.040 s)']);
disp(['Overshoot: ', num2str(info.Overshoot), '% (Req: < 16%)']);
disp(['Steady-State Error: ', num2str(ss_error), '% (Req: 0%)']);
Output :

Task 5:
Code:
%a
m1 = 2500;
m2 = 320;
k1 = 80000;
k2 = 500000;
b1 = 350;
b2 = 15020;
%b
A=[0 1 0 0
-(b1*b2)/(m1*m2) 0 ((b1/m1)*((b1/m1)+(b1/m2)+(b2/m2)))-(k1/m1) -(b1/m1)
b2/m2 0 -((b1/m1)+(b1/m2)+(b2/m2)) 1
k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0];
B=[0 0
1/m1 (b1*b2)/(m1*m2)
0 -(b2/m2)
(1/m1)+(1/m2) -(k2/m2)];
C=[0 0 1 0];
D=[0 0];
sys=ss(A,B,C,D);
%c
Aa = [[A,[0 0 0 0]'];[C, 0]];
Ba = [B;[0 0]];
Ca = [C,0];
Da = D;
sys=ss(Aa,Ba,Ca,Da);
K = [0 2.3e6 5e8 0 8e6]
t = 0:0.01:2;
sys_cl = ss(Aa-Ba(:,1)*K,-0.1*Ba,Ca,Da);
step(sys_cl*[0;1],t)
title('Closed-Loop Response to a 0.1-m Step')
% d: Compute Overshoot and Settling Time
% Simulate response to step input on second input channel (road disturbance)
u = repmat([0; 0.1], 1, length(t)); % 0.1-m step disturbance on input 2
[y, t_out] = lsim(sys_cl, u', t);
% Compute step response characteristics
info = stepinfo(y, t_out);

disp('--- Performance Metrics ---');

disp(['Overshoot: ', num2str(info.Overshoot), ' % (Requirement: < 5%)']);
disp(['Settling Time: ', num2str(info.SettlingTime), ' s (Requirement: < 5s)']);
Output:

Task 6:
Code:
% Inverted Pendulum Control
M = 0.5; % Cart mass
m = 0.2; % Pendulum mass
b = 0.1; % Friction coefficient
I = 0.006; % Pendulum inertia
g = 9.8; % Gravity
l = 0.3; % Pendulum length
den = I*(M+m)+M*m*l^2;
A = [0 1 0 0;
0 -(I+m*l^2)*b/den m^2*g*l^2/den 0;
0 0 0 1;
0 -m*l*b/den m*g*l*(M+m)/den 0];
B = [0; (I+m*l^2)/den; 0; m*l/den];
C = [1 0 0 0; 0 0 1 0];
D = [0; 0];
sys = ss(A,B,C,D);
p = eig(A);
co = ctrb(sys);
cr = rank(co);
Q = C'*C;
R = 1;
K = lqr(A,B,Q,R);
Ac = A-B*K;
sys_cl = ss(Ac,B,C,D);
t = 0:0.01:5;
r = 0.2*ones(size(t));
[y,t,x] = lsim(sys_cl,r,t);
figure;
[ax,h1,h2] = plotyy(t,y(:,1),t,y(:,2),'plot');
ylabel(ax(1),'Position (m)');
ylabel(ax(2),'Angle (rad)');
title('LQR Response');
Q(1,1) = 5000;
Q(3,3) = 100;
K = lqr(A,B,Q,R);
Ac = A-B*K;
sys_cl = ss(Ac,B,C,D);
[y,t,x] = lsim(sys_cl,r,t);
figure;
[ax,h1,h2] = plotyy(t,y(:,1),t,y(:,2),'plot');
ylabel(ax(1),'Position (m)');
ylabel(ax(2),'Angle (rad)');
title('Improved LQR');
Cn = [1 0 0 0];
sys_ss = ss(A,B,Cn,0);
N = rscale(sys_ss,K);
sys_cl = ss(Ac,B*N,C,D);
[y,t,x] = lsim(sys_cl,r,t);
figure;
[ax,h1,h2] = plotyy(t,y(:,1),t,y(:,2),'plot');
ylabel(ax(1),'Position (m)');
ylabel(ax(2),'Angle (rad)');
title('With Precompensation');
ob = obsv(sys_ss);
orb = rank(ob);
cp = eig(Ac);
op = [-40 -41 -42 -43];
L = place(A',Cn',op)';
Ace = [A-B*K B*K; zeros(size(A)) A-L*Cn];
Bce = [B*N; zeros(size(B))];
Cce = [C zeros(size(C))];
Dce = [0;0];
sys_est = ss(Ace,Bce,Cce,Dce);
[y,t,x] = lsim(sys_est,r,t);
figure;
[ax,h1,h2] = plotyy(t,y(:,1),t,y(:,2),'plot');
ylabel(ax(1),'Position (m)');
ylabel(ax(2),'Angle (rad)');
title('With Observer');
function Nbar = rscale(sys,K)
[A,B,C,D] = ssdata(sys);
s = size(A,1);
Z = [A B; C 0];
N = [zeros(s,1); 1];
X = Z\N;
Nx = X(1:s);
Nu = X(end);
Nbar = Nu + K*Nx;
end

Output:
Task 7:
Code:
%% a) Controllability Check
A = [-0.313 56.7 0; -0.0139 -0.426 0; 0 56.7 0];
B = [0.232; 0.0203; 0];
C = [0 0 1];
D = [0];
co = ctrb(A,B);
Controllability = rank(co) % Should be 3 (full rank) for controllable system

%% b) Initial LQR Design

p = 2;
Q = p*C'*C; % Weight on output (pitch angle)
R = 1; % Weight on control input
[K] = lqr(A,B,Q,R);

sys_cl = ss(A - B*K, B, C, D);

step(0.2*sys_cl) % 0.2 radian step input (~11.5 degrees)
ylabel('pitch angle (rad)');
title('Initial LQR Response');

%% c) Improved LQR Design

p = 50; % Increased weight on output
Q = p*C'*C;
R = 1;
[K] = lqr(A,B,Q,R);

sys_cl = ss(A-B*K, B, C, D);

step(0.2*sys_cl)
ylabel('pitch angle (rad)');
title('Improved LQR Response');

%% d) LQR with Precompensation

Nbar = rscale(A,B,C,D,K); % Precompensation gain

sys_cl = ss(A-B*K, B*Nbar, C, D);

step(0.2*sys_cl)
ylabel('pitch angle (rad)');
title('Final LQR Response with Precompensation');
% Required rscale function (add at the end)
function Nbar = rscale(A,B,C,D,K)
s = size(A,1);
Z = [zeros(1,s) 1];
N = inv([A B; C D])*Z';
Nx = N(1:s);
Nu = N(s+1:end);
Nbar = Nu + K*Nx;
end

Output: