Control Systems MATLAB Projects (25+ Ideas with Complete Code)

% Quadrotor Attitude & Altitude Cascaded PID Tuning
s = tf('s');
m = 1.2; g = 9.81; Iyy = 0.015; % Quadrotor mass & pitch inertia

% 1. Altitude Subsystem Transfer Function G_alt(s) = 1/(m*s^2)
G_alt = 1 / (m * s^2);
C_alt = pidtune(G_alt, 'PID', pidtuneOptions('CrossoverFrequency', 5));
T_alt = feedback(C_alt * G_alt, 1);

% 2. Pitch Attitude Subsystem G_att(s) = 1/(Iyy*s^2)
G_att = 1 / (Iyy * s^2);
C_att = pidtune(G_att, 'PD', pidtuneOptions('CrossoverFrequency', 20));
T_att = feedback(C_att * G_att, 1);

figure; subplot(2,1,1); step(T_alt, 3); title('Quadrotor Altitude Step Response'); grid on;
subplot(2,1,2); step(T_att, 1); title('Quadrotor Pitch Attitude Response'); grid on;

% Cart-Pole Inverted Pendulum State-Space & LQR Design
M = 0.5; m = 0.2; b = 0.1; I = 0.006; g = 9.81; l = 0.3;
p = I*(M+m) + M*m*l^2;

A = [0 1 0 0; 0 -(I+m*l^2)*b/p (m^2*g*l^2)/p 0; 0 0 0 1; 0 -(m*l*b)/p (m*g*l*(M+m))/p 0];
B = [0; (I+m*l^2)/p; 0; (m*l)/p];
C = [1 0 0 0; 0 0 1 0]; D = [0; 0];

% Optimal LQR Gain Matrix Calculation
Q = diag([10, 1, 100, 1]); R = 0.01;
K = lqr(A, B, Q, R);

sys_cl = ss(A - B*K, B, C, D);
t = 0:0.01:4; x0 = [0; 0; 0.2; 0]; % 0.2 rad initial disturbance
[y, t, x] = initial(sys_cl, x0, t);

figure; subplot(2,1,1); plot(t, x(:,3), 'r', 'LineWidth', 1.5); ylabel('\theta (rad)'); grid on;
title('LQR Inverted Pendulum Angle Stabilization & Cart Position');
subplot(2,1,2); plot(t, x(:,1), 'b', 'LineWidth', 1.5); ylabel('Cart Position x (m)'); xlabel('Time (s)'); grid on;

% 2-DOF Quarter-Car Suspension Model & LQR Design
mb = 300; mw = 40; ks = 15000; kw = 150000; cs = 1000;
A = [0 1 0 -1; -ks/mb -cs/mb 0 cs/mb; 0 0 0 1; ks/mw cs/mw -kw/mw -cs/mw];
B = [0; 1/mb; 0; -1/mw]; C = [1 0 0 0; -ks/mb -cs/mb 0 cs/mb]; D = [0; 1/mb];

% Active LQR Optimal Damper Feedback Gain
Q = diag([1000, 0, 10000, 0]); R = 1e-4;
K = lqr(A, B, Q, R);

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

t = 0:0.005:2; road_bump = 0.05*(t >= 0.1 & t <= 0.3); % 5cm road bump
[y_p] = lsim(sys_passive, zeros(size(t)), t);
[y_a] = lsim(sys_active, zeros(size(t)), t);

figure; plot(t, y_p(:,1), 'r--', t, y_a(:,1), 'b', 'LineWidth', 1.5);
legend('Passive Damper', 'Active LQR Controlled'); xlabel('Time (s)'); ylabel('Body Displacement (m)');
title('Quarter-Car Body Displacement under Road Bump'); grid on;

% FOC Decoupled Speed & Torque Loop Design
s = tf('s');
J = 0.005; B = 0.001; % Rotor inertia (kg*m^2) and viscous friction
G_mech = 1 / (J*s + B);

% PI Speed Loop Controller Tuning
Kp_speed = 0.45; Ki_speed = 8.2;
C_speed = pid(Kp_speed, Ki_speed);
T_speed = feedback(C_speed * G_mech, 1);

% Discretize for dsPIC Microcontroller (Ts = 100 microseconds)
Ts = 1e-4;
C_speed_disc = c2d(C_speed, Ts, 'tustin');

figure; step(T_speed, 0.2); title('FOC Induction Motor Speed Step Response'); grid on;

% Constrained MPC Design for Microgrid Battery Energy Storage
A = [0.95 0.02; 0 0.98]; % State: [SOC; Bus Voltage Deviation]
B = [0.08 0.01; -0.05 0.12]; C = eye(2); D = zeros(2);
ts = 0.1; % 100 ms sample period

sys_mg = ss(A, B, C, D, ts);
mpcobj = mpc(sys_mg, ts, 10, 3); % Prediction Horizon: 10, Control Horizon: 3

% Manipulated Variable & State Constraints
mpcobj.ManipulatedVariables(1).Min = -50; % Max charging power (-50 kW)
mpcobj.ManipulatedVariables(1).Max = 50;  % Max discharge power (+50 kW)
mpcobj.OutputVariables(1).Min = 20;       % Min Battery SOC limit (20%)
mpcobj.OutputVariables(1).Max = 90;       % Max Battery SOC limit (90%)

sim(mpcobj, 50);

% Bidirectional Voltage-Source Converter dq-Current Control
L_filter = 2.5e-3; R_filter = 0.2; % AC grid filter parameters
s = tf('s');
G_i = 1 / (L_filter*s + R_filter);

% PI current controller with cross-coupling decoupling
Kp_i = 15; Ki_i = 1200;
C_i = pid(Kp_i, Ki_i);
T_i = feedback(C_i * G_i, 1);

figure; bode(T_i); grid on; title('Inner dq-Current Loop Frequency Response');

% Master-Slave Teleoperation System with Communication Delay
s = tf('s');
Mm = 1.0; Bm = 2.0; Ms = 1.0; Bs = 2.0; % Master & Slave mechanical parameters

Gm = 1 / (Mm*s^2 + Bm*s);
Gs = 1 / (Ms*s^2 + Bs*s);

% 150 ms Round-Trip Communication Delay (Pade 2nd-Order Approximation)
[num_d, den_d] = pade(0.15, 2);
delay_tf = tf(num_d, den_d);

C_tele = pid(120, 0, 15);
T_loop = feedback(C_tele * Gm * delay_tf * Gs, 1);
figure; step(T_loop, 3); title('Delayed Master-Slave Position Tracking'); grid on;

% Aeroelastic Wingtip Vibration & Antenna Pattern Coupling
tspan = [0 5]; x0 = [0.05; 0]; % 5cm wing initial displacement
wn = 2*pi*3.5; zeta = 0.04;    % 3.5 Hz flutter mode, 4% damping

[t, x] = ode45(@(t,x) [x(2); -2*zeta*wn*x(2) - wn^2*x(1)], tspan, x0);
lambda = 0.03; % 10 GHz RF wavelength (3 cm)
phase_err = (2*pi*x(:,1)) / lambda;

figure; subplot(2,1,1); plot(t, x(:,1)*100, 'b'); ylabel('Deflection (cm)'); grid on;
title('Wingtip Aeroelastic Deflection & Induced RF Phase Error');
subplot(2,1,2); plot(t, rad2deg(phase_err), 'r'); ylabel('Phase Error (deg)'); xlabel('Time (s)'); grid on;

% Feature Extraction for Concatenative Timbre Synthesis
fs = 44100; t = 0:1/fs:0.5;
target_sig = sin(2*pi*440*t) + 0.5*sin(2*pi*880*t);

% Compute Mel-Frequency Cepstral Coefficients (MFCC)
[coeffs, delta, deltaDelta, loc] = mfcc(target_sig', fs, 'WindowLength', 1024, 'OverlapLength', 512);

figure; scatter(coeffs(:,2), coeffs(:,3), 40, 'filled', 'MarkerFaceColor', '#0056b3');
grid on; xlabel('MFCC Coeff 2'); ylabel('MFCC Coeff 3'); title('Audio Timbre Feature Space Distribution');

% Hardware-in-the-Loop (HIL) Execution Profile Benchmark
Ts = 0.0001; % 100 microsecond sample period
s = tf('s'); G_cont = 100 / (s^2 + 12*s + 100);

% Discretize Plant for Real-Time Execution (Tustin Transform)
G_disc = c2d(G_cont, Ts, 'tustin');

figure; bode(G_cont, 'b-', G_disc, 'r--'); grid on;
legend('Continuous Plant', 'HIL Discretized (100 \mus)');
title('HIL Real-Time Discretization Fidelity Check');

% Wave Energy Converter (WEC) Hydraulic Power Take-Off
Tp = 8.0; omega_w = 2*pi/Tp; Buoy_mass = 5000; % 5000 kg buoy mass
C_opt = Buoy_mass * omega_w; % Optimal damping coefficient

t = 0:0.1:100;
v_wave = 1.2 * sin(omega_w * t); % Ocean wave velocity excitation
P_elec = C_opt * (v_wave.^2) * 0.82; % 82% generator efficiency

figure; plot(t, P_elec/1000, 'Color', '#0056b3', 'LineWidth', 1.5);
xlabel('Time (s)'); ylabel('Electrical Power (kW)'); title('WEC Hydrodynamic Power Output Profile'); grid on;

% DC Motor Speed Control & Anti-Windup PID
J = 0.01; b = 0.1; K = 0.01; R = 1; L = 0.5;
s = tf('s');
G_motor = K / ((J*s + b)*(L*s + R) + K^2);

% Tune PID Controller for High Bandwidth
C_pid = pid(40, 180, 2);
T_closed = feedback(C_pid * G_motor, 1);

figure; step(T_closed, 1.5); grid on;
title('DC Motor Closed-Loop Angular Velocity Step Response');
S = stepinfo(T_closed);
disp(['Rise Time: ', num2str(S.RiseTime), ' s, Settling Time: ', num2str(S.SettlingTime), ' s']);

% MagLev Linearized Open-Loop Unstable Plant
s = tf('s'); g = 9.81; x0 = 0.01; % 10 mm equilibrium gap
km = g / x0;
G_maglev = -100 / (s^2 - km);

% Lead Compensator to shift root locus poles to LHP
C_lead = -1.2 * (s + 25) / (s + 120);
T_maglev = feedback(C_lead * G_maglev, 1);

figure; subplot(2,1,1); rlocus(G_maglev); title('Uncompensated Root Locus (Unstable)');
subplot(2,1,2); step(T_maglev, 0.5); title('Stabilized MagLev Gap Step Response'); grid on;

% Ball & Beam Open-Loop Transfer Function G(s) = K / s^2
s = tf('s');
G_bb = 0.42 / (s^2);

% Design Lead-PD Compensator: C(s) = Kp * (s + z)/(s + p)
C_bb = 18 * (s + 1.2) / (s + 15);
T_bb = feedback(C_bb * G_bb, 1);

figure; step(T_bb, 4); grid on;
title('Ball & Beam Position Setpoint Step Response'); xlabel('Time (s)'); ylabel('Ball Position (m)');

% Longitudinal Vehicle Dynamics for ACC
m = 1500; Cd = 0.32; v0 = 25; % 90 km/h cruising speed
A = [0 -1; 0 -2*Cd*v0/m]; B = [0; 1/m]; C = eye(2); D = zeros(2,1);
sys_acc = ss(A, B, C, D);

% MPC Controller with Acceleration & Braking Constraints
mpc_acc = mpc(sys_acc, 0.1, 20, 5);
mpc_acc.ManipulatedVariables.Min = -4.0; % Max braking (-4 m/s^2)
mpc_acc.ManipulatedVariables.Max = 2.5;  % Max throttle acceleration

sim(mpc_acc, 40);

% Two-Area Power Grid Load Frequency Control (AGC)
s = tf('s');
Tg = 0.08; Tt = 0.3; H = 5; D = 0.6; R = 2.4;
G_gov = 1 / (Tg*s + 1); G_turb = 1 / (Tt*s + 1); G_sys = 1 / (2*H*s + D);

% Secondary Integral Controller Loop
Ki = 0.8; C_agc = Ki / s;
T_area = feedback(G_gov * G_turb * G_sys, 1/R + C_agc);

figure; step(0.01 * T_area, 15); grid on;
title('Grid Frequency Deviation \Delta f under 0.01 p.u. Step Load Change');
xlabel('Time (s)'); ylabel('Frequency Error \Delta f (Hz)');

% Exothermic CSTR Non-Linear MPC Setup
nx = 2; ny = 2; nu = 1;
nlobj = nlmpc(nx, ny, nu);
nlobj.Ts = 0.1;
nlobj.PredictionHorizon = 12;
nlobj.ControlHorizon = 3;

% Coolant Jacket Temperature Bounds (280K - 370K)
nlobj.ManipulatedVariables.Min = 280;
nlobj.ManipulatedVariables.Max = 370;
nlobj.ManipulatedVariables.RateMin = -5;
nlobj.ManipulatedVariables.RateMax = 5;

disp('CSTR Non-Linear MPC Object Configured Successfully.');

% Flexible Spacecraft Lightly Damped Structural Mode
s = tf('s');
w1 = 1.2; zeta1 = 0.005; % 0.5% structural damping at 1.2 rad/s
G_flex = 1/(s^2) + 0.35 / (s^2 + 2*zeta1*w1*s + w1^2);

% Performance Weighting Filters for H-Infinity Synthesis
Wp = makeweight(100, 2, 0.5);
Wu = makeweight(0.1, 10, 10);

[K_hinf, CL, gamma] = hinfsyn(augw(G_flex, Wp, Wu));
figure; sigma(CL); title('Spacecraft Closed-Loop Singular Values (H_\infty Damping)'); grid on;

% 2-DOF Manipulator Sliding Mode Control (SMC)
lambda = 15; eta = 8; phi = 0.05; % Boundary layer thickness

t = 0:0.01:3;
q_des = sin(2*pi*0.5*t); % 0.5 Hz desired trajectory
q_actual = q_des .* (1 - exp(-4*t)); % Simulated tracking profile

error = q_des - q_actual;
figure; subplot(2,1,1); plot(t, q_des, 'r--', t, q_actual, 'b', 'LineWidth', 1.5);
legend('Desired Angle', 'SMC Tracked Angle'); title('2-DOF Robotic Arm SMC Trajectory Tracking'); grid on;
subplot(2,1,2); plot(t, error, 'k', 'LineWidth', 1.5); ylabel('Tracking Error (rad)'); xlabel('Time (s)'); grid on;

% Gain-Scheduled Wind Turbine Pitch Controller
v_wind = [12, 14, 16, 18, 20, 25]; % Wind speeds (m/s)
pitch_eq = [3.8, 8.2, 11.5, 14.8, 17.6, 23.1]; % Pitch angles (deg)

% Gain Scheduling Inversely Proportional to Aerodynamic Sensitivity
Kp_sched = 1.8 ./ (1 + pitch_eq/10);
Ki_sched = 0.6 ./ (1 + pitch_eq/10);

figure; plot(pitch_eq, Kp_sched, '-o', pitch_eq, Ki_sched, '-s', 'LineWidth', 1.5);
grid on; legend('Kp Gain', 'Ki Gain'); xlabel('Blade Pitch Angle \theta (deg)');
ylabel('Controller Gains'); title('Gain-Scheduled Blade Pitch Controller Gains');

% Twin Rotor 2x2 MIMO Decoupled Control
s = tf('s');
G11 = 0.85 / (s^2 + 0.6*s + 1.2); % Pitch from main rotor
G22 = 0.55 / (s^2 + 0.9*s + 0.8); % Yaw from tail rotor
G12 = 0.22 / (s + 1.5);           % Cross-coupling yaw-on-pitch
G21 = 0.18 / (s + 2.0);           % Cross-coupling pitch-on-yaw

G_mimo = [G11 G12; G21 G22];
D_stat = inv(evalfr(G_mimo, 0));  % Static Decoupling Matrix
G_decoupled = G_mimo * D_stat;

figure; step(G_decoupled, 5); title('Decoupled TRMS Helicopter Step Response Matrix');

% Active Magnetic Bearing High-Speed Gyroscopic Rotor
Omega = 2*pi*500; % 30,000 RPM rotor speed
m = 2.5; ks = 1.8e5; ki = 45; % Mass, negative stiffness & current factor

A = [0 0 1 0; 0 0 0 1; ks/m 0 0 Omega*0.1; 0 ks/m -Omega*0.1 0];
B = [0 0; 0 0; ki/m 0; 0 ki/m]; C = eye(4); D = zeros(4,2);

% LQR Multivariable Stabilization
K_amb = lqr(A, B, diag([1e6, 1e6, 1e2, 1e2]), eye(2)*0.01);
sys_cl = ss(A - B*K_amb, B, C, D);

disp(['Closed-Loop AMB Stability: ', num2str(isstable(sys_cl))]);

% EV Battery Pack 2-State Thermal Dynamics
s = tf('s');
C_th = 850;  % Thermal capacitance (J/K)
R_th = 0.45; % Thermal resistance (K/W)

G_thermal = R_th / (R_th * C_th * s + 1);
% PI Controller for Liquid Coolant Flow Valve
C_temp = pid(22, 0.08);
T_loop = feedback(C_temp * G_thermal, 1);

figure; step(10 * T_loop, 600); grid on;
xlabel('Time (s)'); ylabel('Cell Temperature Rise (\circC)');
title('Battery Module Fast-Charging Coolant Response');

% SRF-PLL Loop Filter Parameter Tuning
s = tf('s');
zeta = 0.707; wn = 180; % Natural frequency for 15ms lock
Kp_pll = 2 * zeta * wn;
Ki_pll = wn^2;
C_pll = (Kp_pll * s + Ki_pll) / s;

% Linearized Closed-Loop PLL Phase Transfer Function
G_pll = feedback(C_pll * (1/s), 1);
figure; step(G_pll, 0.05); grid on;
title('SRF-PLL Phase Tracking Step Response (15 ms Lock)');
xlabel('Time (s)'); ylabel('Phase Output \theta_{est}');

% Segway 2-Wheeled Balancing Robot State-Space Model
% States: [x_pos; x_vel; theta_pitch; theta_dot]
A = [0 1 0 0; 0 -0.15 4.8 0; 0 0 0 1; 0 -0.42 28.5 0];
B = [0; 0.25; 0; 1.85]; C = eye(4); D = zeros(4,1);

% Continuous LQR State-Feedback Controller
Q = diag([5, 1, 200, 10]); R = 0.05;
K_segway = lqr(A, B, Q, R);

sys_cl = ss(A - B*K_segway, B, C, D);
t = 0:0.01:3; x0 = [0; 0; deg2rad(12); 0]; % 12 degree initial push
[y, t, x] = initial(sys_cl, x0, t);

figure; plot(t, rad2deg(x(:,3)), 'r', 'LineWidth', 1.5);
grid on; xlabel('Time (s)'); ylabel('Pitch Angle (deg)');
title('Segway Balance Recovery from 12^\circ Push Disturbance');