Power Electronics MATLAB Projects (18+ Ideas with Simulink Models)

% Synchronous Buck Converter Design Parameters
Vin = 48;           % Input DC Voltage (V)
Vout = 12;          % Target Output DC Voltage (V)
fs = 100e3;         % Switching Frequency (100 kHz)
Iout_max = 10;      % Maximum Load Current (10 A)
R_load = Vout / Iout_max; % 1.2 Ohms

% Calculate Minimum Inductance & Capacitance
D = Vout / Vin;     % Nominal Duty Cycle = 0.25
delta_IL = 0.2 * Iout_max; % 20% Inductor Ripple Current
L = (Vin - Vout) * D / (fs * delta_IL); % 45 uH
delta_Vout = 0.01 * Vout;  % 1% Voltage Ripple
C = delta_IL / (8 * fs * delta_Vout);   % 20.8 uF -> Select 47 uF

% Continuous Small-Signal Transfer Function (Control to Output)
Gvd = tf([Vin], [L*C, L/R_load, 1]);

% Tune Closed-Loop PID Controller
C_pid = pidtune(Gvd, 'PID', 2*pi*5000); % 5 kHz Bandwidth
T_cl = feedback(C_pid * Gvd, 1);

figure; step(T_cl); grid on;
title('Closed-Loop Buck Converter Step Response');
fprintf('Designed L: %.2f uH, C: %.2f uF, Gain Margin: %.2f dB\n', L*1e6, C*1e6, allmargin(C_pid*Gvd).GainMargin);

% Perturb and Observe (P&O) MPPT MATLAB Core Function
function D = mppt_po(Vpv, Ipv, delta_D)
    persistent Vprev Pprev Dprev
    if isempty(Vprev)
        Vprev = 0; Pprev = 0; Dprev = 0.5;
    end
    
    Pcurr = Vpv * Ipv;
    delta_P = Pcurr - Pprev;
    delta_V = Vpv - Vprev;
    
    if delta_P > 0
        if delta_V > 0
            D = Dprev - delta_D; % Move towards MPP
        else
            D = Dprev + delta_D;
        end
    elseif delta_P < 0
        if delta_V > 0
            D = Dprev + delta_D;
        else
            D = Dprev - delta_D;
        end
    else
        D = Dprev;
    end
    
    % Saturate Duty Cycle
    D = max(0.1, min(0.9, D));
    Vprev = Vpv; Pprev = Pcurr; Dprev = D;
end

% 3-Phase SPWM Synthesis & Harmonic THD Spectrum
fs = 100e3;  % Sampling Rate (100 kHz)
f0 = 50;     % Fundamental Output Frequency (50 Hz)
fc = 10e3;   % Carrier Switching Frequency (10 kHz)
t = 0:1/fs:0.04; % 2 cycles of 50 Hz
ma = 0.85;   % Modulation Index

% Generate Reference Signals (120 deg shifted)
v_ref_a = ma * sin(2*pi*f0*t);
v_carrier = sawtooth(2*pi*fc*t, 0.5); % Triangular Carrier

% Inverter Switching State (Phase A)
s_a = double(v_ref_a > v_carrier);
Vdc = 400; v_pole_a = Vdc * (s_a - 0.5);

% Compute Total Harmonic Distortion (THD)
thd_inverter = thd(v_pole_a, fs, f0);
fprintf('Total Harmonic Distortion (THD): %.2f dB (%.2f%%)\n', thd_inverter, 10^(thd_inverter/20)*100);

% Space Vector Dwell Time Calculation Function
function [T1, T2, T0, sector] = calc_svpwm(Valpha, Vbeta, Vdc, Ts)
    Vref = sqrt(Valpha^2 + Vbeta^2);
    theta = atan2(Vbeta, Valpha);
    if theta < 0, theta = theta + 2*pi; end
    
    % Sector Identification (1 to 6)
    sector = floor(theta / (pi/3)) + 1;
    theta_rel = theta - (sector - 1) * (pi/3);
    
    % Dwell Time Calculations
    T1 = (sqrt(3) * Ts / Vdc) * Vref * sin(pi/3 - theta_rel);
    T2 = (sqrt(3) * Ts / Vdc) * Vref * sin(theta_rel);
    T0 = Ts - T1 - T2; % Zero Vector Time
end

% Adhesion Vacuum Chamber & BLDC Motor Power Control
m_robot = 4.5;       % Robot Mass (kg)
g = 9.81;            % Gravity
mu = 0.6;            % Friction Coefficient
F_gravity = m_robot * g;
F_adhesion_req = 1.5 * (F_gravity / mu); % Required Adhesion Force

% Vacuum Motor Speed Transfer Function (PWM Duty to Pressure)
s = tf('s');
G_suction = 120 / (0.05*s + 1); % Pressure Pa per % Duty
C_suction = pidtune(G_suction, 'PI');

% Power Consumption Model
V_batt = 24; % 24V LiPo
I_motor = 5.2; % BLDC Suction Current (A)
P_total = V_batt * I_motor;
fprintf('Required Adhesion: %.2f N, Total Suction Power: %.2f W\n', F_adhesion_req, P_total);

% GA-Optimized Fuzzy LFC Objective Function
function itae_score = lfc_fitness(params)
    % params vector encodes FIS membership points
    fis = mamfis('Name', 'LFC_Controller');
    fis = addInput(fis, [-0.5 0.5], 'Name', 'FreqError');
    fis = addMF(fis, 'FreqError', 'trimf', params(1:3), 'Name', 'Negative');
    fis = addMF(fis, 'FreqError', 'trimf', params(4:6), 'Name', 'Zero');
    fis = addMF(fis, 'FreqError', 'trimf', params(7:9), 'Name', 'Positive');
    
    % Simulate Interconnected Power Grid System
    assignin('base', 'opt_fis', fis);
    simOut = sim('two_area_lfc_model', 'StopTime', '20');
    t = simOut.tout;
    e_f = simOut.freq_error;
    
    % Compute Integral Time Absolute Error (ITAE)
    itae_score = trapz(t, t .* abs(e_f));
end

% MATLAB Serial Communication with dsPIC30F4011 Motor Drive
port = "COM3"; baudrate = 115200;
sp = serialport(port, baudrate);
configureTerminator(sp, "LF");

% Send Speed Setpoint Command to dsPIC (e.g. 1200 RPM)
speed_setpoint_rpm = 1200;
writeline(sp, sprintf("SET_RPM:%d", speed_setpoint_rpm));

% Read Live Feedback Stream
flush(sp);
data_line = readline(sp);
motor_data = sscanf(data_line, "V=%f,I=%f,RPM=%f");
fprintf('Measured RPM: %.1f, Current: %.2f A\n', motor_data(3), motor_data(2));

% Symbolic Normal Form Transformation for Generator Swing Equation
syms delta omega M D_d Pm Pmax sym_u
% Non-linear swing dynamics: M*d2(delta)/dt2 + D_d*d(delta)/dt = Pm - Pmax*sin(delta) + u
f_delta = omega;
f_omega = (Pm - Pmax*sin(delta) - D_d*omega + sym_u)/M;

% 3rd Order Taylor Approximation around Operating Point
taylor_omega = taylor(f_omega, [delta, omega], [asin(Pm/Pmax), 0], 'Order', 4);
disp('Higher-Order Modal Decoupling Terms:');
disp(taylor_omega);

% Frequency Decoupled Energy Management Strategy
fc_cutoff = 0.5; % 0.5 Hz Low-Pass Cutoff for Fuel Cell
[b, a] = butter(2, fc_cutoff / (100/2), 'low');

% Demand Power Profile (UDDS Drive Cycle)
t = 0:0.01:100;
P_demand = 30e3 + 15e3 * sin(2*pi*0.1*t) + 10e3 * (t > 30 & t < 40); % Watts

% Filter Base Power to Fuel Cell and Transient to Supercapacitor
P_fuelcell = filter(b, a, P_demand);
P_supercap = P_demand - P_fuelcell; % Handles rapid surge

figure; plot(t, P_demand/1e3, 'k--', t, P_fuelcell/1e3, 'b', t, P_supercap/1e3, 'r');
legend('Total Demand (kW)', 'PEMFC Steady Power', 'Supercap Transient Power'); grid on;

% Accelerometer Step Peak Detection & Stride Integration
fs = 100; % 100 Hz IMU Sampling Rate
t = (0:length(accel_data)-1)/fs;
accel_mag = sqrt(sum(accel_data.^2, 2)) - 9.81; % Remove Gravity

% Detect Footfalls
[pks, step_locs] = findpeaks(accel_mag, 'MinPeakHeight', 1.5, 'MinPeakDistance', 0.4*fs);
step_count = length(step_locs);
stride_length = 0.72; % Meters

fprintf('Detected %d steps. Total distance: %.2f meters.\n', step_count, step_count * stride_length);

% Active PFC Boost Converter Inner Current Loop Design
L_boost = 1e-3;       % 1 mH Boost Inductor
C_dc = 470e-6;        % 470 uF DC Output Cap
V_dc = 400;           % Target DC Bus
Rsense = 0.05;        % Current Shunt Resistor

% Current Loop Plant: Gi(s) = Vdc / (s * L)
s = tf('s');
Gi = (V_dc * Rsense) / (s * L_boost);

% Design PI Current Controller for 15 kHz Crossover
C_current = pidtune(Gi, 'PI', 2*pi*15000);
fprintf('PFC Inner Current Controller Kp: %.4f, Ki: %.2f\n', C_current.Kp, C_current.Ki);

% DAB Theoretical Power Transfer Calculation
V1 = 400; V2 = 800; n = 0.5; % Transformer turns ratio
fs = 100e3; L_leakage = 25e-6; % 25 uH
phi = linspace(-pi/2, pi/2, 200); % Phase Shift Angle (rad)

% SPS Power Transfer Equation
P_dab = (n * V1 * V2) / (2 * pi * fs * L_leakage) .* phi .* (1 - abs(phi)/pi);

figure; plot(phi * 180/pi, P_dab / 1e3, 'LineWidth', 2); grid on;
xlabel('Phase Shift Angle \phi (Degrees)'); ylabel('Transferred Power (kW)');
title('Dual Active Bridge Bidirectional Power vs Phase Shift');

% LLC Resonant Voltage Gain via First Harmonic Approximation
k = 5; % Inductance Ratio Lm/Lr
Q_vals = [0.2, 0.5, 1.0, 2.0]; % Quality factors
fn = linspace(0.5, 2.0, 500);  % Normalized Frequency fsw/fr

figure; hold on;
for Q = Q_vals
    % LLC Gain Magnitude
    M = 1 ./ sqrt( (1 + 1/k - 1./(k*fn.^2)).^2 + Q^2 .* (fn - 1./fn).^2 );
    plot(fn, M, 'LineWidth', 1.8, 'DisplayName', sprintf('Q = %.1f', Q));
end
grid on; xlabel('Normalized Frequency f_n'); ylabel('Voltage Gain M');
title('LLC Resonant Converter Gain Curves'); legend show;

% PMSM Parameters
Rs = 0.45; Ld = 3.2e-3; Lq = 3.2e-3; % Henries
psi_m = 0.12; p = 4; % 4 Pole Pairs

% d-q Current Plant: G(s) = 1 / (L*s + R)
s = tf('s');
G_current = 1 / (Ld*s + Rs);

% Inner Loop PI Controller Tuning for 1 kHz Bandwidth
omega_bw = 2*pi*1000;
Kp_i = Ld * omega_bw;
Ki_i = Rs * omega_bw;
C_pi_current = tf([Kp_i, Ki_i], [1, 0]);

fprintf('Designed FOC Current Loop Kp: %.4f, Ki: %.2f\n', Kp_i, Ki_i);

% 3-Level Phase Disposition (PD) PWM Switching Logic
Vdc = 600; ma = 0.9; f0 = 50; fc = 5000;
t = 0:1/100e3:0.04;
v_ref = ma * sin(2*pi*f0*t);

% Upper and Lower Triangular Carriers in Phase
carrier_upper = sawtooth(2*pi*fc*t, 0.5); % Range 0 to 1
carrier_lower = carrier_upper - 1;       % Range -1 to 0

% 3-Level Pole Voltage Generation
s1 = double(v_ref > carrier_upper);
s2 = double(v_ref > carrier_lower);
v_pole = (Vdc/2) * (s1 + s2 - 1); % Levels: +Vdc/2, 0, -Vdc/2

thd_npc = thd(v_pole, 100e3, f0);
fprintf('3-Level NPC Line Voltage THD: %.2f%%\n', 10^(thd_npc/20)*100);

% SST High-Frequency Isolation Link Sizing
S_rated = 100e3; % 100 kVA Power Rating
f_hf = 20e3;     % 20 kHz High-Frequency Link
B_max = 0.25;    % Nanocrystalline Core Flux Density (Tesla)
Ac = 45e-4;      % Core Cross-sectional Area (m^2)

% Primary Turns Calculation (Faraday's Law)
V_pri = 2000; % 2 kV HF AC Square Wave
N_pri = round(V_pri / (4 * f_hf * B_max * Ac));
fprintf('Required HF Transformer Primary Turns: %d turns (Weight: ~12 kg vs 250 kg 50Hz)\n', N_pri);

% Vienna Rectifier Voltage & Current Loops
V_grid_rms = 230; % 230V Phase RMS
V_dc_total = 800; % 800V (+400V / -400V split bus)
P_out = 15e3;     % 15 kW Telecom Load

% Input Boost Inductor Sizing
fs = 50e3; % 50 kHz Switching Frequency
delta_I = 0.15 * (P_out / (3*V_grid_rms)); % 15% ripple
L_vienna = (V_dc_total/2) / (6 * fs * delta_I);

fprintf('Vienna Rectifier Boost Inductors: %.2f uH each phase\n', L_vienna*1e6);

% SRF d-q Voltage Sag Detection and Injection Control
function [V_inj_a, V_inj_b, V_inj_c] = dvr_controller(Va, Vb, Vc, theta)
    % Clarke Transformation abc -> alpha-beta
    Valpha = (2/3) * (Va - 0.5*Vb - 0.5*Vc);
    Vbeta = (2/3) * (sqrt(3)/2*Vb - sqrt(3)/2*Vc);
    
    % Park Transformation alpha-beta -> dq
    Vd = Valpha*cos(theta) + Vbeta*sin(theta);
    Vq = -Valpha*sin(theta) + Vbeta*cos(theta);
    
    % Reference Target: Vd_ref = 1.0 p.u. (325V peak), Vq_ref = 0
    Vd_ref = 325.2; Vq_ref = 0;
    Vd_inj = Vd_ref - Vd; % Sag Compensating Voltage
    Vq_inj = Vq_ref - Vq;
    
    % Inverse Park Transformation to abc
    V_alpha_inj = Vd_inj*cos(theta) - Vq_inj*sin(theta);
    V_beta_inj = Vd_inj*sin(theta) + Vq_inj*cos(theta);
    
    V_inj_a = V_alpha_inj;
    V_inj_b = -0.5*V_alpha_inj + (sqrt(3)/2)*V_beta_inj;
    V_inj_c = -0.5*V_alpha_inj - (sqrt(3)/2)*V_beta_inj;
end