How to simulate a variable load in an electric microgrids problem?

Hi everyone, I'm tryng to simulate a control problem applied to an electrical microgrid modeled as a network with 2 inverters and 1 load. The equation we used are the following:

I solved the DAE system for a constant load:

clear
clc

%% 1.1. Our equation variables and the parameters are:
syms p_1(t) p_2(t) theta_1(t) theta_2(t) theta_3(t)  P_star_1 P_star_2 P_star_3 D_1 D_2 k_1 k_2 A_12 A_13 A_21 A_23 A_31 A_32 A_23 Lc_12 Lc_21

%% 2.1. The equations we want to represent are:
eqn1 = diff(p_1(t)) == P_star_1/k_1 - p_1(t)/k_1 - (1/k_1)*A_12*sin(theta_1(t)-theta_2(t)) - (1/k_1)*A_13*sin(theta_1(t)-theta_3(t)) - (1/k_1)*Lc_12*(p_1(t)/D_1 - p_2(t)/D_2);
eqn2 = diff(p_2(t)) == P_star_2/k_2 - p_2(t)/k_2 - (1/k_2)*A_21*sin(theta_2(t)-theta_1(t)) - (1/k_2)*A_23*sin(theta_2(t)-theta_3(t)) - (1/k_2)*Lc_21*(p_2(t)/D_2 - p_1(t)/D_1);
eqn3 = diff(theta_1(t)) == P_star_1/D_1 - p_1(t)/D_1 - (1/D_1)*A_12*sin(theta_1(t)-theta_2(t)) - (1/D_1)*A_13*sin(theta_1(t)-theta_3(t));
eqn4 = diff(theta_2(t)) == P_star_2/D_2 - p_2(t)/D_2 - (1/D_2)*A_21*sin(theta_2(t)-theta_1(t)) - (1/D_2)*A_23*sin(theta_2(t)-theta_3(t));
eqn5 = P_star_3-A_31*sin(theta_3(t)-theta_1(t))-A_32*sin(theta_3(t)-theta_2(t))==0;

DAEs = [eqn1; eqn2; eqn3; eqn4; eqn5]

%% 3.1 Place the variables as colomn vector:
DAEvars = [p_1(t); p_2(t); theta_1(t); theta_2(t); theta_3(t)];
origVars = length(DAEvars);

%% 4.1 Check Incidence of Variables:
M_inc = incidenceMatrix(DAEs,DAEvars);

%% 5.1 Verifying DAE index:
isLowIndexDAE(DAEs,DAEvars);

%% 6.1 Defining parameters contained in DAE equations:
pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs,pDAEvars);

% Extra parameters are : A_12, A_13, A_21, A_23, A_31, A_32, D_1, D_2, ...
% Lc_12, Lc_21, P_star_1, P_star_2, P_star_3, k_1, k_2

%% 7. Creating DAE function g:
g = daeFunction(DAEs,DAEvars, A_12, A_13, A_21, A_23, A_31, A_32, D_1, D_2, Lc_12, Lc_21, ...
    P_star_1, P_star_2, P_star_3, k_1, k_2);

%% 8. Parametrization and handling function G:
P_star_1 = 2;
P_star_2 = 3;
%P_star_3 = -4;
P_star_3 = -3;

k_1 = 10e-06;
k_2 = 10e-06;
D_1 = 4000;
D_2 = 6000; 
A_12=1;
A_21=1;
A_13=2.5;
A_31=2.5;
A_32=1.7;
A_23=1.7;
Lc_21 = -1000;
Lc_12 = -1000;

G = @(t,Y,YP) g(t,Y,YP, A_12, A_13, A_21, A_23, A_31, A_32, D_1, D_2, ...
    Lc_12, Lc_21, P_star_1, P_star_2, P_star_3, k_1, k_2);

%% 9. Looking for consistent initial condition:
y0est = [0; 0; 0; 0; 0];
yp0est = zeros(5,1);
[yp, yp0] = decic(G,0, y0est,[], yp0est, []);

%% 10. Risolving DAE with ode15i:
[tSol,ySol] = ode15i(G,[0 6],yp,yp0);

%% The power active injection at each inverteris:
plot(tSol, A_12*sin(ySol(:,3)-ySol(:,4))+A_13*sin(ySol(:,3)-ySol(:,5)),'LineWidth',2) % P_e,1
hold on
plot(tSol, A_21*sin(ySol(:,4)-ySol(:,3))+A_23*sin(ySol(:,4)-ySol(:,5)),'LineWidth',2) % P_e,2
axis([0 6 0.5 3.5])

I don't know how to implement a variable load that change in the interval [2s;4s] because of a changing of the parameter P_star_3 from the value -3 to -4.

The resulting plot must have the following form:

I hope you can help me, thank you for your attention!

[Sam Chak]

@Roberto Palmese, can you sketch for clarity? Do you mean a sudden discontinuous spike?

x = linspace(0, 6, 6001);
y = (sign(x - 2) - sign(x - 4))/2 + 3;
plot(x, y, 'linewidth', 1.5)
grid on
xlabel('t')
ylabel('P_{\ast3}')
ylim([2 5])

[Simone Mormile]

yes, he mean a sudden discontinous spike like that

Answer 1

[Kausthub]

Hi Robert,

I understand that you would like to create a plot with value of “P_star_3” equal to -4 for an interval of [2,4] and with a value of -3 for the rest of the intervals. To achieve this, you could try plotting the graph twice, one with for the interval [2,4] with "P_start_3" equal to -4 and the other with "P_start_3" equals to -3 for the rest of the intervals. To get the values for the [2,4] range you could try conditional indexing. For example.

idx = (tSol<=2) | (tSol>=4);
tSol = tSol(idx);
ySol = ySol(idx,:);

This gives you the values of X and Y for the interval of [2,4] for the plot. Similarly you could plot the second graph with “P_start_3” equal to -3 for the rest of the interval. For your reference, I have plotted the first graph for the [2,4] interval using conditional indexing:

clear
clc
%% 1.1. Our equation variables and the parameters are:
syms p_1(t) p_2(t) theta_1(t) theta_2(t) theta_3(t)  P_star_1 P_star_2 P_star_3 D_1 D_2 k_1 k_2 A_12 A_13 A_21 A_23 A_31 A_32 A_23 Lc_12 Lc_21
%% 2.1. The equations we want to represent are:
eqn1 = diff(p_1(t)) == P_star_1/k_1 - p_1(t)/k_1 - (1/k_1)*A_12*sin(theta_1(t)-theta_2(t)) - (1/k_1)*A_13*sin(theta_1(t)-theta_3(t)) - (1/k_1)*Lc_12*(p_1(t)/D_1 - p_2(t)/D_2);
eqn2 = diff(p_2(t)) == P_star_2/k_2 - p_2(t)/k_2 - (1/k_2)*A_21*sin(theta_2(t)-theta_1(t)) - (1/k_2)*A_23*sin(theta_2(t)-theta_3(t)) - (1/k_2)*Lc_21*(p_2(t)/D_2 - p_1(t)/D_1);
eqn3 = diff(theta_1(t)) == P_star_1/D_1 - p_1(t)/D_1 - (1/D_1)*A_12*sin(theta_1(t)-theta_2(t)) - (1/D_1)*A_13*sin(theta_1(t)-theta_3(t));
eqn4 = diff(theta_2(t)) == P_star_2/D_2 - p_2(t)/D_2 - (1/D_2)*A_21*sin(theta_2(t)-theta_1(t)) - (1/D_2)*A_23*sin(theta_2(t)-theta_3(t));
eqn5 = P_star_3-A_31*sin(theta_3(t)-theta_1(t))-A_32*sin(theta_3(t)-theta_2(t))==0;
DAEs = [eqn1; eqn2; eqn3; eqn4; eqn5]

DAEs = 

%% 3.1 Place the variables as colomn vector:
DAEvars = [p_1(t); p_2(t); theta_1(t); theta_2(t); theta_3(t)];
origVars = length(DAEvars);
%% 4.1 Check Incidence of Variables:
M_inc = incidenceMatrix(DAEs,DAEvars);
%% 5.1 Verifying DAE index:
isLowIndexDAE(DAEs,DAEvars);
%% 6.1 Defining parameters contained in DAE equations:
pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs,pDAEvars);
% Extra parameters are : A_12, A_13, A_21, A_23, A_31, A_32, D_1, D_2, ...
% Lc_12, Lc_21, P_star_1, P_star_2, P_star_3, k_1, k_2
%% 7. Creating DAE function g:
g = daeFunction(DAEs,DAEvars, A_12, A_13, A_21, A_23, A_31, A_32, D_1, D_2, Lc_12, Lc_21, ...
    P_star_1, P_star_2, P_star_3, k_1, k_2);
%% 8. Parametrization and handling function G:
P_star_1 = 2;
P_star_2 = 3;
P_star_3 = -4;
% P_star_3 = -3;
k_1 = 10e-06;
k_2 = 10e-06;
D_1 = 4000;
D_2 = 6000; 
A_12=1;
A_21=1;
A_13=2.5;
A_31=2.5;
A_32=1.7;
A_23=1.7;
Lc_21 = -1000;
Lc_12 = -1000;
G = @(t,Y,YP) g(t,Y,YP, A_12, A_13, A_21, A_23, A_31, A_32, D_1, D_2, ...
    Lc_12, Lc_21, P_star_1, P_star_2, P_star_3, k_1, k_2);
%% 9. Looking for consistent initial condition:
y0est = [0; 0; 0; 0; 0];
yp0est = zeros(5,1);
[yp, yp0] = decic(G,0, y0est,[], yp0est, []);
%% 10. Risolving DAE with ode15i:
[tSol,ySol] = ode15i(G,[0 6],yp,yp0);
%% The power active injection at each inverteris:
idx = (tSol>=2) & (tSol<=4);
tSol = tSol(idx);
ySol = ySol(idx,:);
% newySol = 
plot(tSol, A_12*sin(ySol(:,3)-ySol(:,4))+A_13*sin(ySol(:,3)-ySol(:,5)),'LineWidth',2) % P_e,1
hold on
plot(tSol, A_21*sin(ySol(:,4)-ySol(:,3))+A_23*sin(ySol(:,4)-ySol(:,5)),'LineWidth',2) % P_e,2
axis([0 6 0.5 3.5])

Hope this helps!