Matlab: изменение параметра ODE онлайн с условием IF

Question

Итак, я пытаюсь решить 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, уменьшалось на 5. Так что-то вроде:

if Vsol23s << 0.95*E
a = a - 5;
end

Я знаючто сейчас Ode сначала решает для Gsol23s, а затем вычисляет Vsol23s.Есть ли способ обойти?

Любые идеи высоко ценятся.

---

РЕДАКТИРОВАТЬ 1: Решение, которое дает один шаг за один раз (а не вся матрица для tspan),мог бы добиться цели, но я не знаю, возможно ли это.

Answer 1

Параметризация вашего difq затем передать обновленные значения 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

Answer 2

Ну, как я уже говорил, я не думаю, что есть способ сделать пошаговое решение, используя функцию ODExy Matlab. Я решаю проблему с помощью RK4, который имеет почти такое же решение, как ODExy, но, поскольку он не адаптивный, есть гораздо больше точек. Я представляю код, чтобы его могли использовать и другие.

Чтобы прояснить проблему, у нас есть система электроснабжения с индуктивной линией передачи без потерь (X = 10 Ом, R = 0 Ом) и нагрузкой. Параллельно этой нагрузке есть несколько конденсаторов, которые можно использовать, замыкая или открывая соответствующие переключатели. Таким образом, если мощность, потребляемая нагрузкой, достаточно высока, чтобы снизить напряжение нагрузки (VL) ниже 95% от напряжения источника (VS), заглушка подключается для увеличения VL путем компенсации реактивной мощности. Если VL> 105% VS, крышка отключается.

Д.Е. что дает поведение нагрузки: 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');

Надеюсь, это поможет!