0

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.

  • Use events and the solution continuation method. Keep the ODE function simple and smooth. – LutzL Apr 16 at 8:14
1

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
  • This "kind of" works, but want to hold the values of Vsol23s and Psol23s when Vsol23s < 0.95*E. Eg. the first 3 values (with a = 0) give Vsol23s > 0.95*E. The forth value (which makes Vsol23s < 0.95*E) must be calculated for a = -5, and with initial conditions the values of the third step. I am starting to believe that there is no way to do this with odeXX matlab function and that I need to program my own arithmetic method (possibly RK4). Anyway, thanks for your time! – thece Apr 16 at 16:33
  • @thece I edited the above to possibly reflect what you're looking for - here I stop the ode solution, update a, and restart the solution at new values of t0 & G0. – peterfranciscook Apr 16 at 16:59
0

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!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.