Matlab: Changing ODE Parameter Online with IF-Condition

0
votes

So, I am trying to solve the ODE that describes the response and stability of a power system. The precise equation and theory does not really matter. As the title suggests, I want the ode solver to change some parameter according to an IF condition. Here is the code:

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');

I want, every time that, V is less than 0.95*E, to decrease a by 5. So something like:

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

I know that right now, ode first solves for Gsol23s and then calculates the Vsol23s. Is there any way around?

Any ideas are much appreciated.

EDIT 1: A solution that gives one ode step at a time (and not the whole matrix for tspan), could do the trick, but I don't know if this is possible.

2
matlabode
Use events and the solution continuation method. Keep the ODE function simple and smooth. - Lutz Lehmann

2 Answers

1
votes

Parameterize your diffeq then pass updated values of a to it as needed. 

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
0
votes

Well, as I said before, I don't think that there is a way to do a step-by-step solution, using ODExy Matlab function. I solve the problem using RK4, which has almost exactly the same solution as ODExy, but since its not adaptive, there are way more points. I present the code, so it can be used by others too.

To clarify the problem, we have a electrical power system, with a inductive lossless transmission line (X = 10 Ω, R = 0 Ω) and a load. Parallel to that load there are multiple capacitors that can be used by closing or opening the corresponding switches. So, if the power drawn by the load is high enough to drop the load voltage (VL) under 95% of source voltage (VS), a cap is connected in order to increase the VL, by compensating the reactive power. If VL > 105% VS, the cap is disconnected.

The D.E. that gives the load's behavior is 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');

Hope it helps!