Matlab：使用 IF 条件在线更改 ODE 参数答案

【问题标题】：Matlab: Changing ODE Parameter Online with IF-ConditionMatlab：使用 IF 条件在线更改 ODE 参数
【发布时间】：2019-04-24 17:50: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 的整个矩阵）的解决方案可以解决问题，但我不知道这是否可能。

【问题讨论】：

使用事件和解决方案延续方法。保持 ODE 函数简单流畅。

标签： matlab ode

【解决方案1】：

参数化您的 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;
end

【讨论】：

这种“类型”有效，但希望在 Vsol23s 0.95*E。第四个值（使 Vsol23s
@thece 我编辑了上面的内容以可能反映您正在寻找的内容 - 在这里我停止 ode 解决方案，更新 a，然后以 t0 和 G0 的新值重新启动解决方案。

【解决方案2】：

好吧，正如我之前所说，我不认为有一种方法可以使用 ODExy Matlab 函数进行逐步解决。我使用 RK4 解决问题，它的解决方案与 ODExy 几乎完全相同，但由于它不是自适应的，因此有更多的点。我把代码展示出来，其他人也可以使用。

为了澄清这个问题，我们有一个电力系统，带有一个感应无损传输线（X = 10 Ω，R = 0 Ω）和一个负载。与该负载并联，有多个电容器可以通过关闭或打开相应的开关来使用。因此，如果负载汲取的功率高到足以使负载电压 (VL) 下降到电源电压 (VS) 的 95% 以下，则连接一个电容器以通过补偿无功功率来增加 VL。如果 VL > 105% VS，则封盖断开。

D.E.给出负载的行为是 dG/dt = P0/T - G(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');

希望对你有帮助！

【讨论】：