clear; clc; close all;
%% ========= 目标/工况 =========
Vin = 100; % 100Vdc
Vdc_ref = 400; % 400Vdc
P_rated = 700; % 700W
f_grid = 50; % 50Hz
w_grid = 2*pi*f_grid;
% 三相线电压 220Vrms
Vll_rms = 220;
Vph_rms = Vll_rms/sqrt(3);
Vph_pk = sqrt(2)*Vph_rms;
%% ========= 电路参数（平均模型） =========
plant = struct();
% LCL（每相）
plant.L1 = 1.5e-3; plant.R1 = 0.10; % inverter-side
plant.L2 = 0.8e-3; plant.R2 = 0.10; % grid-side
plant.Cf = 4.7e-6; plant.Rd = 1.0; % series Rd + Cf branch (阻尼)
% Boost + DC link
plant.Lb = 1.0e-3;
plant.RLb = 0.05;
plant.Cdc = 1000e-6;
plant.Rbleed = 5000;
%% ========= 仿真设置 =========
dt = 5e-6; % 电路积分步长
Ts_ctrl = 50e-6; % 20kHz 控制更新
Nctrl = max(1, round(Ts_ctrl/dt));
t_end = 0.40;
N = floor(t_end/dt) + 1;
t = (0:N-1)'*dt;
x = zeros(11,1);
x(1) = 0.0; % IL
x(2) = 360.0; % Vdc 初始
%% ========= 日志 =========
log.Vdc = zeros(N,1);
log.IL = zeros(N,1);
log.mabc = zeros(N,3);
log.dboost= zeros(N,1);
log.vgabc = zeros(N,3);
log.igabc = zeros(N,3);
log.theta = zeros(N,1);
log.omega = zeros(N,1);
log.idref = zeros(N,1);
log.iqref = zeros(N,1);
log.vd = zeros(N,1);
log.vq = zeros(N,1);
%% ========= 控制量保持 =========
u = struct();
u.d = 0.0;
u.m = [0;0;0];
theta = 0; omega = w_grid; id_ref = 0; iq_ref = 0; vd = 0; vq = 0;
%% ========= 主循环 =========
for k = 1:N
tk = t(k);
% 三相电网相电压（相对中性点）
va = Vph_pk*sin(w_grid*tk);
vb = Vph_pk*sin(w_grid*tk - 2*pi/3);
vc = Vph_pk*sin(w_grid*tk + 2*pi/3);
vg = [va; vb; vc];
% 并网电流（取 L2 后的 i2，方向=逆变器->电网）
ia = x(9); ib = x(10); ic = x(11);
% 功率给定：0.05s 后上电到 700W
if tk < 0.05
Pcmd = 0;
else
Pcmd = P_rated;
end
% 每 Ts_ctrl 更新一次"单片机算法"
if mod(k-1, Nctrl) == 0
[d_boost, m_a, m_b, m_c, theta, omega, id_ref, iq_ref, vd, vq] = ...
mcu_control_3ph_step( ...
x(2), x(1), Vin, va, vb, vc, ia, ib, ic, Pcmd);
m0 = (m_a + m_b + m_c)/3;
m_a = m_a - m0; m_b = m_b - m0; m_c = m_c - m0;
u.d = d_boost;
u.m = [m_a; m_b; m_c];
end
% 记录
log.Vdc(k) = x(2);
log.IL(k) = x(1);
log.mabc(k,:) = u.m.';
log.dboost(k) = u.d;
log.vgabc(k,:) = vg.';
log.igabc(k,:) = [ia ib ic];
log.theta(k) = theta;
log.omega(k) = omega;
log.idref(k) = id_ref;
log.iqref(k) = iq_ref;
log.vd(k) = vd;
log.vq(k) = vq;
if k == N, break; end
% RK4 积分
x = rk4_step(@(xx) plant_f(xx, u, vg, Vin, plant), x, dt);
end
%% ========= 计算 THD & PF（末尾10个基波周期） =========
res = calc_thd_pf(t, log.vgabc, log.igabc, f_grid);
fprintf('\n===== 结果（取末尾10个基波周期）=====\n');
fprintf('PF_total = %.4f\n', res.pf_total);
fprintf('THD_A = %.2f %%\n', res.thd_a);
fprintf('THD_B = %.2f %%\n', res.thd_b);
fprintf('THD_C = %.2f %%\n', res.thd_c);
fprintf('P_total = %.1f W\n', res.P_total);
fprintf('Q_total = %.1f var\n', res.Q_total);
fprintf('S_total = %.1f VA\n', res.S_total);
%% ========= 画图 =========
figure; plot(t, log.Vdc, 'LineWidth', 1.2); grid on;
xlabel('t (s)'); ylabel('Vdc (V)'); title('DC Link Voltage (should ~400V)');
figure; plot(t, log.dboost, 'LineWidth', 1.2); grid on;
xlabel('t (s)'); ylabel('d_{boost}'); title('Boost Duty');
figure;
plot(t, log.igabc(:,1), 'LineWidth', 1.0); hold on;
plot(t, log.igabc(:,2), 'LineWidth', 1.0);
plot(t, log.igabc(:,3), 'LineWidth', 1.0);
grid on; xlabel('t (s)'); ylabel('i_{grid} (A)'); title('Grid Currents i_a,i_b,i_c');
legend('ia','ib','ic');
figure;
plot(t, log.vgabc(:,1), 'LineWidth', 1.0); hold on;
plot(t, log.vgabc(:,2), 'LineWidth', 1.0);
plot(t, log.vgabc(:,3), 'LineWidth', 1.0);
grid on; xlabel('t (s)'); ylabel('v_{grid} (V)'); title('Grid Voltages v_a,v_b,v_c');
legend('va','vb','vc');
figure;
plot(t, log.vq, 'LineWidth', 1.2); grid on;
xlabel('t (s)'); ylabel('v_q (V)'); title('PLL alignment indicator (vq -> 0)');
%% ===================== 子函数区 =====================
function xnext = rk4_step(f, x, h)
k1 = f(x);
k2 = f(x + 0.5*h*k1);
k3 = f(x + 0.5*h*k2);
k4 = f(x + h*k3);
xnext = x + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
end
function dx = plant_f(x, u, vg, Vin, p)
% 平均电路模型：Boost + DC Link + 三相两电平逆变 + LCL(等效Rd串Cf) + 电网源
IL = x(1);
Vdc = max(x(2), 10);
i1 = x(3:5); % inverter-side inductor currents
vC = x(6:8); % capacitor voltage
i2 = x(9:11); % grid-side inductor currents (= grid current, 方向=逆变器->电网)
d = min(max(u.d, 0.0), 0.95);
m = u.m;
% 逆变器相电压（平均模型）：v_inv = 0.5*m*Vdc
v_inv = 0.5 .* m .* Vdc;
% LCL 结点：i_f = i1 - i2
i_f = i1 - i2;
% 结点电压 v_node = vC + Rd*i_f
v_node = vC + p.Rd .* i_f;
% L1: L1 di1/dt = v_inv - v_node - R1*i1
di1 = (v_inv - v_node - p.R1.*i1) ./ p.L1;
% L2: L2 di2/dt = v_node - v_grid - R2*i2
di2 = (v_node - vg - p.R2.*i2) ./ p.L2;
% Cf: dvC/dt = i_f / C
dvC = i_f ./ p.Cf;
% Boost: Lb dIL/dt = Vin - (1-d)*Vdc - RLb*IL
dIL = (Vin - (1-d)*Vdc - p.RLb*IL) / p.Lb;
% DC link: Cdc dVdc/dt = i_boost_to_dc - i_dc_inv - Vdc/Rbleed
i_boost_to_dc = (1-d)*IL;
% 逆变器从直流侧取电流：i_dc = p_ac / Vdc（p_ac = sum(v_inv * i1)）
p_inv = sum(v_inv .* i1); % 逆变器AC侧瞬时功率（>0 表示逆变器向AC输出）
i_dc_inv = p_inv / Vdc;
dVdc = (i_boost_to_dc - i_dc_inv - Vdc/p.Rbleed) / p.Cdc;
dx = [dIL; dVdc; di1(:); dvC(:); di2(:)];
end
%% ======== "单片机逻辑"控制========
function [d_boost, m_a, m_b, m_c, theta, omega, id_ref, iq_ref, vd, vq] = ...
mcu_control_3ph_step(Vdc, IL, Vin, va, vb, vc, ia, ib, ic, P_cmd)
Ts_ctrl = 50e-6; % 20kHz
Ts_outer = 500e-6; % 2kHz
persistent p st
if isempty(p)
p = default_params();
st = default_state();
end
% 1) Boost 双环（Vdc->ILref->duty）
[st, d_boost] = boost_control_step(st, p, Vdc, IL, Vin, Ts_ctrl, Ts_outer);
% 2) SRF-PLL（对三相电压）
[st, theta, omega, vd, vq] = srf_pll_step(st, p, va, vb, vc, Ts_ctrl);
% 3) 功率->电流给定（电压对齐：P=(3/2)*vd*id，iq=0）
vd_eff = max(abs(vd), p.vd_min);
id_ref = -(2.0 * P_cmd) / (3.0 * vd_eff);
iq_ref = 0.0;
id_ref = sat(id_ref, -p.i_grid_max, p.i_grid_max);
% 4) 电流采样->dq
[i_alpha, i_beta] = clarke(ia, ib, ic);
[id_meas, iq_meas] = park(i_alpha, i_beta, theta);
% 5) dq PI + 解耦 + 电压前馈 -> v_alpha_beta_ref
[st, v_alpha_ref, v_beta_ref] = current_dq_control_step( ...
st, p, id_ref, iq_ref, id_meas, iq_meas, vd, vq, omega, Ts_ctrl);
% 6) SVPWM -> duty -> modulation(-1..1)
[da, db, dc] = svpwm_duty(v_alpha_ref, v_beta_ref, Vdc);
m_a = 2.0 * da - 1.0;
m_b = 2.0 * db - 1.0;
m_c = 2.0 * dc - 1.0;
end
function p = default_params()
p = struct();
p.Vdc_ref = 400.0;
p.w_nom = 2*pi*50;
p.vd_min = 20.0;
p.i_grid_max = 8.0;
% 控制器用等效电感/电阻（按 L1+L2 量级）
p.Leq = 2.5e-3;
p.R_eq = 0.2;
% PLL
p.pll_kp = 60.0;
p.pll_ki = 4000.0;
p.pll_omega_min = 2*pi*45;
p.pll_omega_max = 2*pi*55;
% 电流环（带宽 ~300Hz 起步）
wc_i = 2*pi*300;
p.i_kp = p.Leq * wc_i;
p.i_ki = p.R_eq * wc_i;
p.vref_max_scale = 0.95;
% Boost 外环（Vdc->ILref）
p.boost_Imax = 20.0;
p.boost_Imin = 0.0;
p.boost_outer_kp = 0.08;
p.boost_outer_ki = 40.0;
% Boost 内环（IL->duty）
p.boost_inner_kp = 0.02;
p.boost_inner_ki = 200.0;
p.boost_d_max = 0.95;
end
function st = default_state()
st = struct();
st.vdc_int = 0.0;
st.il_int = 0.0;
st.il_ref = 0.0;
st.outer_div = uint32(0);
st.pll_int = 0.0;
st.theta = 0.0;
st.omega = 2*pi*50;
st.id_int = 0.0;
st.iq_int = 0.0;
end
function [st, d_boost] = boost_control_step(st, p, Vdc, IL, Vin, Ts_ctrl, Ts_outer)
% 外环调度（慢环）
if st.outer_div == 0
ev = p.Vdc_ref - Vdc;
st.vdc_int = sat(st.vdc_int + ev * Ts_outer, -50.0, 50.0);
il_ref = p.boost_outer_kp * ev + p.boost_outer_ki * st.vdc_int;
st.il_ref = sat(il_ref, p.boost_Imin, p.boost_Imax);
end
N = uint32(max(1, round(Ts_outer / Ts_ctrl)));
st.outer_div = st.outer_div + 1;
if st.outer_div >= N
st.outer_div = uint32(0);
end
% 内环（快环）
ei = st.il_ref - IL;
st.il_int = sat(st.il_int + ei * Ts_ctrl, -2.0, 2.0);
% 前馈必须用"当前Vdc"，否则母线会越抬越高
d_ff = 1.0 - Vin / max(Vdc, 1.0);
d_pi = p.boost_inner_kp * ei + p.boost_inner_ki * st.il_int;
d_boost = d_ff + d_pi;
% 当母线高于参考且IL_ref已到0，允许彻底停止升压
if (st.il_ref <= (p.boost_Imin + 1e-6)) && (Vdc > p.Vdc_ref + 2.0)
d_boost = 0.0;
end
d_boost = sat(d_boost, 0.0, p.boost_d_max);
end
function [st, theta, omega, vd, vq] = srf_pll_step(st, p, va, vb, vc, Ts)
[v_alpha, v_beta] = clarke(va, vb, vc);
[vd, vq] = park(v_alpha, v_beta, st.theta);
e = -vq; % 目标 vq -> 0
st.pll_int = sat(st.pll_int + e*Ts, -200.0, 200.0);
omega = p.w_nom + p.pll_kp*e + p.pll_ki*st.pll_int;
omega = sat(omega, p.pll_omega_min, p.pll_omega_max);
theta = wrap_0_2pi(st.theta + omega*Ts);
st.theta = theta;
st.omega = omega;
end
function [st, v_alpha_ref, v_beta_ref] = current_dq_control_step( ...
st, p, id_ref, iq_ref, id_meas, iq_meas, vd, vq, omega, Ts)
ed = id_ref - id_meas;
eq = iq_ref - iq_meas;
st.id_int = sat(st.id_int + ed*Ts, -200.0, 200.0);
st.iq_int = sat(st.iq_int + eq*Ts, -200.0, 200.0);
ud = p.i_kp*ed + p.i_ki*st.id_int;
uq = p.i_kp*eq + p.i_ki*st.iq_int;
% 解耦 + 前馈（Leq 模型）
v_d_ref = vd + ud + p.R_eq*id_meas - omega*p.Leq*iq_meas;
v_q_ref = vq + uq + p.R_eq*iq_meas + omega*p.Leq*id_meas;
[v_alpha_ref, v_beta_ref] = inv_park(v_d_ref, v_q_ref, st.theta);
end
function [da, db, dc] = svpwm_duty(v_alpha, v_beta, Vdc)
Vdc = max(Vdc, 10.0);
Vref = hypot(v_alpha, v_beta);
Vref_max = 0.95 * Vdc / sqrt(3);
if Vref > Vref_max
s = Vref_max / max(Vref, 1e-9);
v_alpha = v_alpha*s; v_beta = v_beta*s;
Vref = Vref_max;
end
ang = atan2(v_beta, v_alpha);
if ang < 0, ang = ang + 2*pi; end
sector = floor(ang/(pi/3)) + 1;
sector = min(max(sector,1),6);
ang_s = ang - (sector-1)*(pi/3);
k = sqrt(3)/Vdc;
T1 = k * Vref * sin(pi/3 - ang_s);
T2 = k * Vref * sin(ang_s);
T0 = 1 - T1 - T2;
if T0 < 0, T0 = 0; end
T0_2 = 0.5*T0;
switch sector
case 1
da = T1 + T2 + T0_2; db = T2 + T0_2; dc = T0_2;
case 2
da = T1 + T0_2; db = T1 + T2 + T0_2; dc = T0_2;
case 3
da = T0_2; db = T1 + T2 + T0_2; dc = T2 + T0_2;
case 4
da = T0_2; db = T1 + T0_2; dc = T1 + T2 + T0_2;
case 5
da = T2 + T0_2; db = T0_2; dc = T1 + T2 + T0_2;
otherwise
da = T1 + T2 + T0_2; db = T0_2; dc = T1 + T0_2;
end
da = sat(da,0,1); db = sat(db,0,1); dc = sat(dc,0,1);
end
function [alpha, beta] = clarke(a,b,c)
alpha = (2/3)*(a - 0.5*b - 0.5*c);
beta = (2/3)*((sqrt(3)/2)*b - (sqrt(3)/2)*c);
end
function [d,q] = park(alpha,beta,theta)
ct = cos(theta); st = sin(theta);
d = alpha*ct + beta*st;
q = -alpha*st + beta*ct;
end
function [alpha,beta] = inv_park(d,q,theta)
ct = cos(theta); st = sin(theta);
alpha = d*ct - q*st;
beta = d*st + q*ct;
end
function y = sat(x,lo,hi)
y = min(max(x,lo),hi);
end
function th = wrap_0_2pi(th)
twopi = 2*pi;
th = th - floor(th/twopi)*twopi;
end
%% ======== THD / PF 计算 ========
function results = calc_thd_pf(t, vabc, iabc, f0)
if nargin < 4, f0 = 50; end
T0 = 1/f0;
t_end = t(end);
t_start = max(t(1), t_end - 10*T0);
idx = (t >= t_start);
tt = t(idx);
vv = vabc(idx,:);
ii = iabc(idx,:);
vv = vv - mean(vv,1);
ii = ii - mean(ii,1);
fs = 1/mean(diff(tt));
N = length(tt);
V = fft(vv).*2/N;
I = fft(ii).*2/N;
f = (0:N-1)'*fs/N;
[~, k1] = min(abs(f - f0));
max_h = 40;
bins = k1*(2:max_h);
bins = bins(bins < N/2);
thd = zeros(1,3);
for ph = 1:3
I1 = abs(I(k1,ph));
Ih = abs(I(bins,ph));
thd(ph) = 100 * sqrt(sum(Ih.^2)) / max(I1,1e-9);
end
Vrms = rms(vv,1);
Irms = rms(ii,1);
p_inst = sum(vv.*ii, 2);
P_total = mean(p_inst);
S_total = sum(Vrms .* Irms);
pf_total = P_total / max(S_total, 1e-9);
% 粗略无功：用基波相量夹角估
V1 = squeeze(V(k1,:)); I1 = squeeze(I(k1,:));
phi = angle(V1) - angle(I1);
Q_total = sum( (Vrms .* Irms) .* sin(phi) );
results = struct();
results.thd_a = thd(1); results.thd_b = thd(2); results.thd_c = thd(3);
results.P_total = P_total;
results.Q_total = Q_total;
results.S_total = S_total;
results.pf_total = pf_total;
end
逐行注释这段代码，并且说明运行结果中的5副图片各代表着什么
订阅专栏 解锁全文
扫一扫
576
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
