matlab代码记录帖

这个帖子主要是给我自己看的，因为不想做重复性的工作，所以代码就记录在这里好了。

下面是弱电网下并网逆变器的稳定性分析及增强控制_unlocked这篇论文里的推导

s=tf('s');
Gi=(2*80*4*pi*s)/(s^2+2*4*pi*s+(100*pi)^2)+2;
Gd=(1-0.75*1/(9.6*1000)*s)/((1+0.75*1/(9.6*1000)*s));
GL=1/(10e-3+(0.25e-3)*s);
Gf=1/((1/2000*2*pi)*s^2+(1/2000*2*pi*0.707)*s+1);
Z=(1+Gi*Gd*GL)/(GL*(Gf*Gd-1))
Zg=0.5e-3*s;
sys=Zg/Z;
nyquist(sys)

%这个是直接简洁的化简代入的，但是好像和paper上还是相差很大
s=tf('s');
Gi=(2009.6*s)/(s^2+25.12*s+314^2)+2;
Gd=(1-(0.078125e-3)*s)/(1+(0.078125e-3)*s);
GL=1/(10e-3+(0.25e-3)*s);
Gf=1/((0.079618e-3)*s^2+0.707*(0.079618e-3)*s+1)
Z=(1+Gi*Gd*GL)/(GL*(Gf*Gd-1))
Zg=0.5e-3*s;
sys=Zg/Z;
nyquist(sys)

%这是paper上给出的
%SCR=11.7;
num=[2.416e-24,8.586e-20,1.529e-5,1.686e-11,1.218e-7,5.741e-4,0.5667,86.57,3474];
den=[9.663e-21,4.201e-16,9.074e-12,681.209e-10,0.073e-4,6.81e-3,0.46,0]
Z=tf(num,den);%Z的传递函数
s=tf('s');
Zg=(0.5e-3)*s;%Zg的传递函数
sys=Zg/Z;
nyquist(sys)

下面是帮师兄计算变压器参数的:

%L1表示原边漏感
%e1是副边开路，e2是副边短路，e3是是原边开路，e4是原边短路
%首先是计算T1电感
syms L1 L2 L3;
N=0.5/1.2;
e1=L1+L3-5.153;
%e2=L1+L2*N^2*L3/(L3+N^2*L2)-0.102;
e3=L2+L3/(N^2)-29.093;
e4=L2+((L1/N^2)*(L3/N^2))/((L1/N^2)+(L3/N^2))-1.046;
[x0,y0,z0] = solve(e1,e3,e4,L1,L2,L3);
%[x1,y1,z1] = solve(e1,e2,e4,L1,L2,L3);
vpa(x0,6),vpa(y0,6),vpa(z0,6)%保留六位有效数字
%vpa(x1,6),vpa(y1,6),vpa(z1,6)

%计算T1电阻
syms L1 L2 L3;
N=0.5/1.2;
e1=L1+L3-215.8;
% e2=L1+L2*N^2*L3/(L3+N^2*L2)-1.85;
e3=L2+L3/(N^2)-1.23*1000;
e4=L2+((L1/N^2)*(L3/N^2))/((L1/N^2)+(L3/N^2))-145.27;
[x0,y0,z0] = solve(e1,e3,e4,L1,L2,L3);%T1计算电阻
% [x1,y1,z1] = solve(e1,e2,e4,L1,L2,L3);
vpa(x0,6),vpa(y0,6),vpa(z0,6)%保留六位有效数字
% vpa(x1,6),vpa(y1,6),vpa(z1,6)

%计算T2电感
syms L1 L2 L3;
N=0.5/1.2;
e1=L1+L3-5.170;
% e2=L1+L2*N^2*L3/(L3+N^2*L2)-0.130;
e3=L2+L3/(N^2)-28.694;
e4=L2+((L1/N^2)*(L3/N^2))/((L1/N^2)+(L3/N^2))-1.08;
[x0,y0,z0] = solve(e1,e3,e4,L1,L2,L3);%T2计算电感
% [x1,y1,z1] = solve(e1,e2,e4,L1,L2,L3);
vpa(x0,6),vpa(y0,6),vpa(z0,6)%保留六位有效数字
% vpa(x1,6),vpa(y1,6),vpa(z1,6)

%计算T2电阻
syms L1 L2 L3;
N=0.5/1.2;
e1=L1+L3-221.83;
% e2=L1+L2*N^2*L3/(L3+N^2*L2)-3.41;
e3=L2+L3/(N^2)-1.2268*1000;
e4=L2+((L1/N^2)*(L3/N^2))/((L1/N^2)+(L3/N^2))-165.33;
[x0,y0,z0] = solve(e1,e3,e4,L1,L2,L3);%T2计算电阻
% [x1,y1,z1] = solve(e1,e2,e4,L1,L2,L3);
vpa(x0,6),vpa(y0,6),vpa(z0,6)%保留六位有效数字
% vpa(x1,6),vpa(y1,6),vpa(z1,6)

LCL并网并网逆变器

%LCL谐振，L1=150uH,L2=50uH,C=100uF
s=tf('s');
num1=[1];
den1=[(150e-6)*(50e-6)*(100e-6),0,(150e-6)+(50e-6),0]%L1=150uH,L2=50uH;C=100uF
sys1=tf(num1,den1);
bode(sys1),hold on

num2=[1];
den2=[200e-6,0]% 单L滤波器
sys2=tf(num2,den2);
bode(sys2),hold on

num3=[1];
den3=[(150e-6)*(50e-6)*(100e-6)*0.1,(200e-6)*(100e-6)*0.1,(150e-6)+(50e-6),0]%电容串电阻，无源阻尼，r=0.1Ω
sys3=tf(num3,den3);
bode(sys3),hold on

%考虑电网阻抗的时候，
figure
s=tf('s');
L1=4e-3;
L2=1e-3;
C=9e-6;
Lg=1e-3;
Rg=0;
for n=1:10:21
G0=1/(s^3*L1*C*(L2+n*Lg)+s^2*L1*C*n*Rg+s*(L2+n*Lg+L1)+n*Rg);
bode(G0),hold on
end