Matlab:用中频条件在线改变ODE参数
Matlab:用中频条件在线改变ODE参数
提问于 2019-04-15 22:54:53
因此，我试图解决描述电力系统响应和稳定性的ODE。精确的方程式和理论并不重要。如标题所示，我希望ode求解器根据IF条件更改一些参数。以下是代码：
clear;
E = 10;
X = 10;
R = 0;
T = 5; % Time Constant
tmax = 100;
G0 = 0; % Initial Value
Psc = E^2/X;
P0 = 1.25*Psc; % Power Losses
G = 0:0.001:2;
a = tand (0);
V = E ./ sqrt( (1+R*G+a*X*G).^2 + (X*G-a*R*G).^2 );
P = E^2*G ./ ( (1+R*G+a*X*G).^2 + (X*G-a*R*G).^2 );
% ODE Numerical Solution
[t23, Gsol23s] = ode23s(@(t, G) P0/T - G* (E ./ sqrt( (1+R*G+a*X*G).^2 ...
+ (X*G-a*R*G)).^2 )^2, [0 tmax], G0);
Vsol23s = E ./ sqrt( (1+R*Gsol23s+a*X*Gsol23s).^2 ...
+ (X*Gsol23s-a*R*Gsol23s).^2 );
Psol23s = E^2*Gsol23s ./ ( (1+R*Gsol23s+a*X*Gsol23s).^2 ...
+ (X*Gsol23s-a*R*Gsol23s).^2 );
figure;
plot(P/Psc, V/E, 'linewidth', 3);
hold on;
grid on;
plot(Psol23s/Psc, Vsol23s/E, 'linewidth', 2);
xlabel('P/Psc'); ylabel('V/E')
figure;
plot(t23, Vsol23s/E, t23, Psol23s/Psc, t23, Gsol23s*X)
legend('V/E', 'P/Psc', 'G*B');我希望，每次，V小于0.95*E时，把a减少5。
if Vsol23s << 0.95*E
a = a - 5;
end我现在知道，ode首先为Gsol23s求解，然后计算Vsol23s。有办法绕过吗？
任何想法都是非常感谢的。
编辑1:一次给出一个ode步骤(而不是tspan的整个矩阵)的解决方案可以做到这一点，但我不知道这是否可能。
回答 2
Stack Overflow用户
发布于 2019-04-16 15:02:37
参数化您的diffeq，然后根据需要将a的更新值传递给它。
E = 10;
X = 10;
R = 0;
tmax = 100;
G0 = 0; % Initial Value
Psc = E^2/X;
P0 = 1.25*Psc; % Power Losses
G = 0:0.001:2;
a = tand (0);
V = E ./ sqrt( (1+R*G+a*X*G).^2 + (X*G-a*R*G).^2 );
P = E^2*G ./ ( (1+R*G+a*X*G).^2 + (X*G-a*R*G).^2 );
[t23, Gsol23s] = ode23s(@(t, G) f(t, G, a), [0 tmax], G0);
Vsol23s = E ./ sqrt( (1+R*Gsol23s+a*X*Gsol23s).^2 + (X*Gsol23s-a*R*Gsol23s).^2 );
Psol23s = E^2*Gsol23s ./ ( (1+R*Gsol23s+a*X*Gsol23s).^2 + (X*Gsol23s-a*R*Gsol23s).^2 );
tIdx = find(Vsol23s < 0.95*E, 1, 'first');
if tIdx
t0 = t23(tIdx);
G0 = Gsol23s(tIdx);
tSol = t23(1:tIdx-1);
gSol = Gsol23s(1:tIdx-1);
vSol = Vsol23s(1:tIdx-1);
pSol = Psol23s(1:tIdx-1);
end
while tIdx
a = a - 5;
[t23, Gsol23s] = ode23s(@(t, G) f(t, G, a), [t0 tmax], G0);
Vsol23s = E ./ sqrt( (1+R*Gsol23s+a*X*Gsol23s).^2 + (X*Gsol23s-a*R*Gsol23s).^2 );
Psol23s = E^2*Gsol23s ./ ( (1+R*Gsol23s+a*X*Gsol23s).^2 + (X*Gsol23s-a*R*Gsol23s).^2 );
tIdx = find(Vsol23s < 0.95*E, 1, 'first');
disp(tIdx)
if tIdx
t0 = t23(tIdx);
G0 = Gsol23s(tIdx);
tSol = [tSol; t23(1:tIdx-1)];
gSol = [gSol; Gsol23s(1:tIdx-1)];
vSol = [vSol; Vsol23s(1:tIdx-1)];
pSol = [pSol; Psol23s(1:tIdx-1)];
else
tSol = [tSol; t23];
gSol = [gSol; Gsol23s];
vSol = [vSol; Vsol23s];
pSol = [pSol; Psol23s];
end
end
figure;
plot(P/Psc, V/E, 'linewidth', 3);
hold on;
grid on;
plot(pSol/Psc, vSol/E, 'linewidth', 2);
xlabel('P/Psc'); ylabel('V/E')
figure;
plot(tSol, vSol/E, tSol, pSol/Psc, tSol, gSol*X)
legend('V/E', 'P/Psc', 'G*B');
function y = f(t, G, a)
% parameters
E = 10;
X = 10;
R = 0;
T = 5; % Time Constant
Psc = E^2/X;
P0 = 1.25*Psc; % Power Losses
y = P0/T - G* (E ./ sqrt( (1+R*G+a*X*G).^2 + (X*G-a*R*G)).^2 )^2;
endStack Overflow用户
发布于 2019-04-17 09:04:08
嗯，就像我前面说过的，我不认为有一种方法一步一步地解决问题，使用ODExy Matlab函数。我使用RK4来解决这个问题，它的解决方案与ODExy几乎完全相同，但是由于它不是自适应的，所以有更多的点。我给出了代码，这样其他人也可以使用它。
为了澄清这个问题，我们有一个电力系统，一个感性无损传输线(X = 10Ω，R=0Ω)和一个负载。与该负载平行的还有多个电容器，可以通过关闭或打开相应的开关来使用。因此，如果负载所产生的功率足够高，使负载电压( VL )降至源电压(VS)的95%以下，则通过补偿无功功率，连接一个上限以增加VL。如果VL > 105% VS，则帽被断开。
给出载荷行为的dG/dt =P0/ The (T)*V^2/T。
clear;
E = 10; % Source Voltage
X = 10; % Line Impedance
Psc = E^2/X; % Short-Circuit Power
P0 = 0.4*Psc; % Power Losses
G = 0:0.001:2; % Load Conductivity
a = 0; % atand(0), load angle
% % Thevenin Theorem before Load Terminals
beta = 0;
Et = E/(1-beta); Xt = X/(1-beta);
betastep = 0.1;
betamat = 0:betastep:4*betastep;
for k = 1 : length(betamat)
Etmat(k) = E/(1-betamat(k));
Xtmat(k) = X/(1-betamat(k));
for j = 1:length(G)
V(k,j) = Etmat(k) ./ sqrt( (1+a*Xtmat(k)*G(j)).^2 ...
+ (Xtmat(k)*G(j)).^2 );
P(k,j) = V(k,j)^2*G(j);
end
end
% % % %
% Runge Kutta Method
% 4th Order
% % % %
tol = 1e-8; tmax = 10;
e = 1; i = 1; h = 0.0001;
T = 5; % Time Constant
Gsol(1) = 0;
t(1) = 0;
Vsol(1) = Et;
Psol(1) = 0;
% % ODE to be solved
dG = @(t,G) (P0 - G*Et^2/((1+a*G*Xt)^2 + (Xt*G)^2));
while e > tol && max(t) < tmax
% % In every cycle, if beta has changed, change the
% % corresponding values
Et = 1/(1-beta)*E; Xt = X/(1-beta);
dG = @(t,G) (P0 - G*Et^2/((1+a*G*Xt)^2 + (Xt*G)^2));
k1 = h*dG(t(i),Gsol(i));
k2 = h*dG(t(i)+h/2, Gsol(i)+k1/2);
k3 = h*dG(t(i)+h/2, Gsol(i)+k2/2);
k4 = h*dG(t(i)+h, Gsol(i)+k3);
Gsol(i+1) = Gsol(i) + 1/6 * (k1 + 2*k2 + 2*k3 + k4);
t (i+1) = t(i) + h;
Vsol(i+1) = Et / sqrt( (1+a*Xt*Gsol(i+1))^2 + (Xt*Gsol(i+1))^2 );
Psol(i+1) = Vsol(i+1)^2 * Gsol(i+1);
% % With Restriction
if Vsol(i+1) <= 0.95*E && beta < .1 && beta >= 0
beta = beta + betastep;
elseif Vsol(i+1) >= 1.05*E && beta < .1 && beta >= 0
beta = beta - betastep;
end
if beta < 0
beta = 0;
end
if beta < 0
beta = 0;
end
e = abs(Gsol(i+1) - Gsol(i));
i = i+1;
end
Vload = linspace(0, max(Vsol), length(Gsol));
figure;
hold on;
for k = 1:length(betamat)
plot(P(k,:)/Psc, V(k,:)/E, 'linewidth', 3);
end
grid on;
plot(Psol/Psc, Vsol/E, 'linewidth', 1.5);
plot(Gsol.*(Vload.^2)/Psc, Vload/E);
plot([P0/Psc P0/Psc], [0 1], 'k')
xlabel('P/Psc'); ylabel('V/E')
figure;
plot(t, Vsol/E, t, Psol/P0, t, Gsol*X)
grid on; legend('V/E', 'P/P0', 'G*X');希望能帮上忙！
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/55698274
复制相关文章
相似问题
腾讯云开发者
