含有双向储能的电力系统小干扰稳定的matlab计算程序

电力系统的稳定性问题一直是电力工程研究的重点之一。含有双向储能的电力系统能够有效地解决电网的功率平衡问题，提高系统的稳定性。在matlab中，可以通过模拟电力系统的动态响应来研究系统的稳定性。下面是一个简单的含有双向储能的电力系统的matlab计算程序，可以用来进行系统稳定性分析。

% 定义系统参数
Sbase = 100; % 基准容量
Vbase = 10; % 基准电压
Zbase = Vbase^2/Sbase; % 基准阻抗
fbase = 50; % 基准频率
wbase = 2*pi*fbase; % 基准角频率

% 定义系统拓扑结构
B = [2 -1 0; -1 3 -1; 0 -1 2]; % 节点导纳矩阵
G = real(B); % 节点电导矩阵
B = imag(B); % 节点电纳矩阵
C = [1 0 0; 0 1 0; 0 0 1]; % 测量矩阵
D = [0 0 0; 0 0 0; 0 0 0]; % 直通矩阵

% 定义负荷和发电机
Pl = [1; 1; 1]*0.5; % 负荷有功功率
Ql = [1; 1; 1]*0.2; % 负荷无功功率
Pg = [1; 1; 1]*1.2; % 发电机有功功率

% 定义储能装置
Pb_max = [1; 1; 1]*0.5; % 储能装置最大输出功率
Pb_min = [1; 1; 1]*(-0.5); % 储能装置最大吸收功率
E_max = [1; 1; 1]*1; % 储能装置最大储能容量
E_min = [1; 1; 1]*(-1); % 储能装置最小储能容量
Pb0 = [1; 1; 1]*0; % 储能装置初始功率
E0 = [1; 1; 1]*0; % 储能装置初始储能容量
eff = [1; 1; 1]*0.95; % 储能装置充放电效率
Pb = Pb0; % 储能装置实时输出功率
E = E0; % 储能装置实时储能容量

% 定义仿真参数
t_start = 0; % 仿真起始时间
t_stop = 10; % 仿真终止时间
t_step = 0.01; % 仿真步长
t = t_start:t_step:t_stop; % 时间矢量
N = length(t); % 仿真步数

% 定义系统状态变量
theta = zeros(3, N); % 节点相角
omega = ones(3, N); % 节点频率
V = ones(3, N)*Vbase; % 节点电压
Pinj = zeros(3, N); % 节点注入有功功率
Qinj = zeros(3, N); % 节点注入无功功率

% 初始状态
theta(:, 1) = [0; -pi/2; pi/2];
omega(:, 1) = [1; 1; 1];
V(:, 1) = [1; 1; 1]*Vbase;

% 仿真
for i = 1:N-1
    % 计算节点注入功率
    Pinj(:, i) = Pg - Pl - Pb(:, i);
    Qinj(:, i) = Qg - Ql;
    
    % 计算节点电流
    I = Y*V(:, i);
    
    % 计算储能装置功率
    Pb(:, i+1) = Pb(:, i) + t_step*(eff.*(Pinj(:, i) - I(1:3).*conj(V(:, i))));
    Pb(:, i+1) = min(max(Pb(:, i+1), Pb_min), Pb_max);
    E(:, i+1) = E(:, i) + t_step*(Pb(:, i) - eff.*Pb(:, i+1));
    E(:, i+1) = min(max(E(:, i+1), E_min), E_max);
    
    % 计算节点相角和频率
    dtheta = omega(:, i)*t_step;
    domega = t_step*(Pinj(:, i) - I(1:3).*conj(V(:, i)) - Pb(:, i));
    theta(:, i+1) = theta(:, i) + dtheta;
    omega(:, i+1) = omega(:, i) + domega/wbase;
    
    % 计算节点电压
    V(:, i+1) = V(:, i).*exp(1i*theta(:, i+1));
end

% 画图
figure;
subplot(2, 1, 1);
plot(t, theta(1, :), 'r', t, theta(2, :), 'g', t, theta(3, :), 'b');
xlabel('Time (s)');
ylabel('Phase Angle (rad)');
title('Node Phase Angle');
legend('Node 1', 'Node 2', 'Node 3');
subplot(2, 1, 2);
plot(t, omega(1, :), 'r', t, omega(2, :), 'g', t, omega(3, :), 'b');
xlabel('Time (s)');
ylabel('Frequency (Hz)');
title('Node Frequency');
legend('Node 1', 'Node 2', 'Node 3');

上述程序中，首先定义了系统的基本参数，包括基准容量、基准电压、基准阻抗、基准频率和基准角频率。然后定义了系统的拓扑结构，包括节点导纳矩阵、测量矩阵和直通矩阵。接着定义了负荷和发电机的有功和无功功率，以及储能装置的最大输出功率、最大吸收功率、最大储能容量、最小储能容量、初始功率、初始储能容量和充放电效率。然后定义了仿真参数，包括仿真起始时间、仿真终止时间、仿真步长和时间矢量。接着定义了系统状态变量，包括节点相角、节点频率、节点电压、节点注入有功功率和节点注入无功功率。然后定义了初始状态。最后进行仿真，并画出节点相角和频率随时间的变化图。

需要注意的是，这个程序只是一个简单的示例程序，实际的电力系统可能会更加复杂，需要根据具体情况进行调整。此外，由于系统中含有储能装置，所以需要特别注意储能装置的充放电过程。