基于matlab实现的电力系统稳定性分析摆幅曲线代码

完整程序：

 clear;
clc;

t = 0;
tf = 0;
tfl = 0.5;
tc = 0.5; % tc = 0.05, 0.125, 0.5 sec for 2.5 cycles, 6.25 cycles & 25 cycles resp
ts = 0.05;
m = 2.52 / (180 * 50);
i = 2;
dt = 21.64 * pi / 180;
ddt = 0;
time(1) = 0;
ang(1) = 21.64;
pm = 0.9;
pm1 = 2.44;
pm2 = 0.88;
pm3 = 2.00;
speed(1) = 0;
power(1) = 0;
voltage(1) = 1.0;

while t < tfl
    switch t
        case tf
            pam = pm - pm1 * sin(dt);
            pap = pm - pm2 * sin(dt);
            paav = (pam + pap) / 2;
            pa = paav;
        case tc
            pam = pm - pm2 * sin(dt);
            pap = pm - pm3 * sin(dt);
            paav = (pam + pap) / 2;
            pa = paav;
        otherwise
            if t > tf && t < tc
                pa = pm - pm2 * sin(dt);
            elseif t > tc
                pa = pm - pm3 * sin(dt);
            end
    end
    
    ddt = ddt + (ts * ts * pa / m);
    dt = (dt * 180 / pi + ddt) * pi / 180;
    dtdg = dt * 180 / pi;
    t = t + ts;
    
    % Calculate additional parameters
    rotor_speed = ddt * 180 / pi;
    elec_power = voltage(i-1) * sin(dt);
    gen_voltage = voltage(i-1) - elec_power * 0.05;
    
    % Store values in arrays for plotting
    time(i) = t;
    ang(i) = dtdg;
    speed(i) = rotor_speed;
    power(i) = elec_power;
    voltage(i) = gen_voltage;
    
    i = i + 1;
end

figure;
subplot(2, 2, 1);
plot(time, ang, 'k+-');
title('转角度');
xlabel('Time (s)');
ylabel('转角度 (degrees)');

subplot(2, 2, 2);
plot(time, speed, 'r*-');
title('转角度');
xlabel('时间(s)');
ylabel('转速(degrees/s)');

subplot(2, 2, 3);
plot(time, power, 'bo-');
title('电力');
xlabel('时间(s)');
ylabel('能量 (pu)');

subplot(2, 2, 4);
plot(time, voltage, 'g^-');
title('发电机端电压');
xlabel('时间 (s)');
ylabel('电压 (pu)');