【数值分析】发电机方程的数值求解及MATLAB仿真

2023-11-29 15:18:46

1 发电机方程

有发电机的二阶微分方程如下：

其中δ为发电机转子角（rotor angle），ω为角速度（rotor angular speed），H为发电机惯性常数（inertia constant），D为发电机阻尼系数（damping constant）。

2 欧拉法

将微分方程离散化，在离散点处用差商代替导数。其中，根据不同的差商，欧拉法分为向前欧拉和向后欧拉。

2.1 向前欧拉法

2.2 向后欧拉法

向后欧拉是一种隐式格式，计算每一个迭代步时都需要求解一个非线性方程。

2.3 两点欧拉公式

两点欧拉公式需要用到两个初值，其中y0给定，y1一般通过向前欧拉算得。

3 欧拉预估-校正法

由于梯形公式是一种隐式格式，实际计算时不是很方便，因此将其修正成显示格式。

4 发电机方程的数值求解

向前欧拉法：

iter = 100;
h = 1/100;

x0 = [x0_af(1);x0_af(4)];
x(:,1) = x0;
omega = x(1,:);
delta = x(2,:);
Y1 = [ 1.0608 - 2.6305i   0.1082 + 0.6064i   0.1572 + 0.9459i];
K=[];

for i = 1:iter

    k1 = f1(omega(i),delta(i),i,Y1);
    k2 = f1(omega(i)+0.5*h*k1(1),delta(i)+0.5*h*k1(2),i,Y1);
    k3 = f1(omega(i)+0.5*h*k2(1),delta(i)+0.5*h*k2(2),i,Y1);
    k4 = f1(omega(i)+h*k3(1),delta(i)+h*k3(2),i,Y1);
    k22 = f1(omega(i)+h*k1(1),delta(i)+h*k1(2),i,Y1);
    %x(:,i+1) = x(:,i) + (h/6)*(k1+2*k2+2*k3+k4); %四阶龙格库塔方法
    x(:,i+1) = x(:,i) + h*k1; %欧拉方法
    % x(:,i+1) = x(:,i) + 0.5*h*(k1+k22); %改进的欧拉方法
    omega = x(1,:);
    delta = x(2,:);
    K = [K,k1];
end

子函数

%% 4-order Runge-Kutta method 求k1 k2 k3 k4的函数
% 在power_grid ode45_RK4中被调用
%%
function k = f1(omega,delta,i,Y1)
k = zeros(2,1);
load Xaf_out
delta1 = [delta,Xaf_out(i+1,5),Xaf_out(i+1,6)];
E_pf =[1.0566,1.0502,1.0170];
clear sum
Pe1 = E_pf(1).*sum(E_pf.*(real(Y1).*cos(delta*ones(1,3)-delta1)+imag(Y1).*sin(delta*ones(1,3)-delta1)));
k(1) = (pi*60/23.6400)*(0.6611 - Pe1);
k(2) = omega - 2*pi*60;
end

5 ode45

ode45是MATLAB自带的解微分方程的函数，用的是四阶龙格库塔方法

[Taf_out , Xaf_out] = ode45(@(t,x)Gen_Fun(t,x,GEN,Yraf,E_pf),tspan,x0_af,options);

子函数：

function Xdf_out = Gen_Fun(t,x,GEN,Ynn,E_pf)
Xdf_out = zeros(6,1);
omega_r = 2*pi*60;
delta_ij_matrix = zeros(3);
for i = 1:3
    for j = 1:3
        delta_ij_matrix(i,j) = x(3+i) - x(3+j);
    end
end
H = GEN(:,4);
Pm = GEN(:,5);
Pe = E_pf*E_pf'.*(real(Ynn).*cos(delta_ij_matrix) + imag(Ynn).*sin(delta_ij_matrix));
Pe = sum(Pe,2);
Xdf_out(1:3) = omega_r*(Pm - Pe)./(2*H);
Xdf_out(4:6) = x(1:3) - omega_r;
end

如上图是利用欧拉法求解，和ode45求解的对比。

码农公寓