%      SyncPowerSyst.M
%   
%  Synchronisation of interconnected power systems by using the theory of
%  Kuramoto oscillators
%    
%  
% J. Lunze
% 
% 17.9.2018
% Version v. 18.9.2018
% for 2nd edition: 18.5.2021
%
echo off
clear
close all

%%    Parameters and models
%
N1=4;   % number of generators
M1=5;   % inertia in s
M2=5;
M3=5;
M4=5;
D1=15*2*pi;  % in MW/Hz = 2pi MWs/rad  15
D2=18*2*pi;                            18
D3=25*2*pi;
D4=14*2*pi;
V1=110;      % in kV
V2=110;
V3=110;
V4=110;
X12=20;         % in Ohm
X23=15;
X24=36;
X34=20;
%   Network
AG=[0 1 0 0;
   1 0 1 1;
   0 1 0 1;
   0 1 1 0];
k12=V1*V2/X12;
k23=V2*V3/X23;
k24=V2*V4/X24;
k34=V3*V4/X34;
AG=[0   k12  0   0;
   k12  0  k23 k24;
   0   k23  0  k34;
   0   k24  k34 0];
%
Ts=0.001;
Tend=10;
Time=[0:Ts:Tend];
kend=length(Time)-1;
%
%  State-space model of the generators
nTS=3;
nGS=N1*nTS;
A1=[0    1   0;
    0 -D1/M1 1;
    0    0   0];
B1=[0; 1/M1; 0];
C1=[1 0 0];
A2=[0    1   0;
    0 -D2/M2 1;
    0    0   0];
B2=[0; 1/M2; 0];
C2=[1 0 0];
A3=[0    1   0;
    0 -D3/M3 1;
    0    0   0];
B3=[0; 1/M3; 0];
C3=[1 0 0];
A4=[0    1   0;
    0 -D4/M4 1;
    0    0   0];
B4=[0; 1/M4; 0];
C4=[1 0 0];
%   discrete-time models
Generator1=ss(A1, B1, C1, 0);
Generator1d=c2d(Generator1, Ts);
A1d=Generator1d.A;
B1d=Generator1d.B;
Generator2=ss(A2, B2, C2, 0);
Generator2d=c2d(Generator2, Ts);
A2d=Generator2d.A;
B2d=Generator2d.B;
Generator3=ss(A3, B3, C3, 0);
Generator3d=c2d(Generator3, Ts);
A3d=Generator3d.A;
B3d=Generator3d.B;
Generator4=ss(A4, B4, C4, 0);
Generator4d=c2d(Generator4, Ts);
A4d=Generator4d.A;
B4d=Generator4d.B;
%  Overall system
AGS=zeros(nGS,nGS);
AGS(1:nTS,1:nTS)=A1d;
AGS(nTS+1:2*nTS,nTS+1:2*nTS)=A2d;
AGS(2*nTS+1:3*nTS,2*nTS+1:3*nTS)=A3d;
AGS(3*nTS+1:4*nTS,3*nTS+1:4*nTS)=A4d;
BGS=zeros(nGS,N1);
BGS(1:nTS,1)=B1d;
BGS(nTS+1:2*nTS,2)=B2d;
BGS(2*nTS+1:3*nTS,3)=B3d;
BGS(3*nTS+1:4*nTS,4)=B4d;
CGS=zeros(N1, nGS);
CGS(1,1:nTS)=C1;
CGS(2,nTS+1:2*nTS)=C2;
CGS(3,2*nTS+1:3*nTS)=C3;
CGS(4,3*nTS+1:4*nTS)=C4;




%%   Behaviour for balanced power generation
%
figure(1)
deltaP1=0;
deltaP2=0;
deltaP3=0;
deltaP4=0;
theta1=0;
theta2=0;
theta3=0;
theta4=0;
%   initial states
theta10=45*pi/180;    % 45 degrees
theta20=20*pi/180;
theta30=30*pi/180;
theta40=4*30*pi/180-(theta10+theta20+theta30);
x10=[theta10; 0; deltaP1/M1];
x20=[theta20; 0; deltaP2/M2];
x30=[theta30; 0; deltaP3/M3];
x40=[theta40; 0; deltaP4/M4];
X0=[x10; x20; x30; x40];   
%
X=[];                           % state as row vector
U=[];
YGS=[];
X(1,:)=X0;
Y=CGS*X0;                       % vector of current phases
YGS(1,:)=Y';                    % output as row vector
for k1=1:kend
    U=zeros(N1,1);
    for k2=1:N1
        for k3=1:N1
            U(k2,1)=U(k2,1)-AG(k2,k3)*sin(Y(k2)-Y(k3));
        end
    end
    X(k1+1,:)=AGS*X(k1,:)'+BGS*U;
    UGS(k1+1,:)=U';
    Y=CGS*X(k1+1,:)';
    YGS(k1+1,:)=Y';
end
subplot(2,1,1)
plot(Time,YGS(:,1:4),'b');
axis([0 1 0.3 0.8]);
latexylabel('$$\tilde{\theta}_i$$ in rad $$\;$$');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 0.2 0.4 0.6 0.8 1]);
yticks([-0.2 0 0.2 0.4 0.6 0.8]);
subplot(2,1,2)
plot(Time,110*sin(YGS(:,1:2)+2*pi*50*Time'),'b')
axis([0 0.15 -120 120]);
latexylabel('$$v_1, v_2$$ in kV');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 0.05 0.1 0.15]);
yticks([-110 0 110]);



%%   Disturbance by a load change
%
figure(2)
deltaP1=0;
deltaP2=0;
deltaP3=0;
deltaP4=0;
%  Disturbance by load change
pL=zeros(kend+1,4);
t1=0.2;     % begin of disturbance  0.2
t2=0.5;     % 0.5
pLMag=50; % 10 MW
for k1=ceil(t1/Ts):ceil(t2/Ts)
    pL(k1,1)=pLMag;
end
%   initial states
theta10=30*pi/180;   
theta20=30*pi/180; 
theta30=30*pi/180; 
theta40=30*pi/180; 
x10=[theta10; 0; deltaP1/M1];
x20=[theta20; 0; deltaP2/M2];
x30=[theta30; 0; deltaP3/M3];
x40=[theta40; 0; deltaP4/M4];
X0=[x10; x20; x30; x40];   
%
X=[];                           % state as row vector
U=[];
YGS=[];
X(1,:)=X0;
Y=CGS*X0;                       % vector of current phases
YGS(1,:)=Y';                    % output as row vector
for k1=1:kend
    U=zeros(N1,1);
    for k2=1:N1
        for k3=1:N1
            U(k2,1)=U(k2,1)-AG(k2,k3)*sin(Y(k2)-Y(k3));
        end
    end
    X(k1+1,:)=AGS*X(k1,:)'+BGS*U+BGS*pL(k1,:)';
    UGS(k1+1,:)=U';
    Y=CGS*X(k1+1,:)';
    YGS(k1+1,:)=Y';
end
subplot(2,1,1)
plot(Time,YGS(:,1:4),'b');
axis([0 1 0.4 0.7]);
latexylabel('$$\tilde{\theta}_i$$ in rad');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 0.2 0.4 0.6 0.8 1]);
yticks([-0.2 0 0.2 0.4 0.6 0.8]);
subplot(2,1,2)
k1=ceil(0.1/Ts);
k2=ceil(0.2/Ts);
TimeV=Time(1,k1:k2)+0.65;   % shifted time
plot(Time, V1*sin(YGS(:,1:2)+2*pi*50*Time'),'b', TimeV, V1*sin(2*pi*50*TimeV'), 'b--');
axis([0.75 0.85 -120 120]);
latexylabel('$$v_1, v_2$$ in kV');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0.75 0.8 0.85]);
yticks([-110 0 110]);


%%   Behaviour for stationary tie line power deltaPi
%
figure(3)
deltaP1=60;   % in MW
deltaP23=69;
deltaP24=100;
deltaP2=deltaP23+deltaP24-deltaP1;  
deltaP34=130;
deltaP3=deltaP34-deltaP23;
deltaP4=-(deltaP1+deltaP2+deltaP3);
%  stationary phases
theta1=45*2*pi/360;    % fix one of the phases, in rad
theta2=theta1-asin(deltaP1/k12);
theta3=theta2-asin(deltaP23/k23);
theta4=theta3-asin(deltaP34/k34);
%  test of stationarity
deltaP1-k12*sin(theta1-theta2)
deltaP2-(k12*sin(theta2-theta1)+k23*sin(theta2-theta3)+k24*sin(theta2-theta4))
deltaP3-(k23*sin(theta3-theta2)+k34*sin(theta3-theta4))
deltaP4-(k24*sin(theta4-theta2)+k34*sin(theta4-theta3))
%
%  Disturbance by load change
pL=zeros(kend+1,4);
t1=0.2;     % begin of disturbance  0.2
t2=0.5;     % 0.5
pLMag=50; % 10 MW
for k1=ceil(t1/Ts):ceil(t2/Ts)
    pL(k1,1)=pLMag;
end
%   initial states
theta10=theta1;   
theta20=theta2;
theta30=theta3;
theta40=theta4;
x10=[theta10; 0; deltaP1/M1];
x20=[theta20; 0; deltaP2/M2];
x30=[theta30; 0; deltaP3/M3];
x40=[theta40; 0; deltaP4/M4];
X0=[x10; x20; x30; x40];   
%
X=[];                           % state as row vector
U=[];
YGS=[];
X(1,:)=X0;
Y=CGS*X0;                       % vector of current phases
YGS(1,:)=Y';                    % output as row vector
for k1=1:kend
    U=zeros(N1,1);
    for k2=1:N1
        for k3=1:N1
            U(k2,1)=U(k2,1)-AG(k2,k3)*sin(Y(k2)-Y(k3));
        end
    end
    X(k1+1,:)=AGS*X(k1,:)'+BGS*U+BGS*pL(k1,:)';
    UGS(k1+1,:)=U';
    Y=CGS*X(k1+1,:)';
    YGS(k1+1,:)=Y';
end
%
subplot(2,1,1)
plot(Time,YGS(:,1:4),'b');
axis([0 1 0.3 0.9]);
latexylabel('$$\tilde{\theta}_i$$ in rad$$\;$$');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 0.2 0.4 0.6 0.8 1]);
yticks([-0.2 0 0.2 0.4 0.6 0.8]);
%
subplot(2,1,2)
k1=ceil(0.1/Ts);
k2=ceil(0.2/Ts);
TimeV=Time(1,k1:k2)+0.65;   % shifted time
plot(Time,V1*sin(YGS(:,1:2)+2*pi*50*Time'),'b', TimeV, V1*sin(2*pi*50*TimeV'), 'b--');
axis([0.75 0.78 -120 120]);
latexylabel('$$v_1, v_2$$ in kV');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0.75 0.76 0.77 0.78]);
yticks([-110 0 110]);


%%   Behaviour for smaller inertia
%
figure(4)
M1=1;   % inertia in s
M2=1;
M3=1;
M4=1;
%  State-space model of the generators
A1=[0    1   0;
    0 -D1/M1 1;
    0    0   0];
A2=[0    1   0;
    0 -D2/M2 1;
    0    0   0];
A3=[0    1   0;
    0 -D3/M3 1;
    0    0   0];
A4=[0    1   0;
    0 -D4/M4 1;
    0    0   0];
%   discrete-time models
Generator1=ss(A1, B1, C1, 0);
Generator1d=c2d(Generator1, Ts);
A1d=Generator1d.A;
Generator2=ss(A2, B2, C2, 0);
Generator2d=c2d(Generator2, Ts);
A2d=Generator2d.A;
Generator3=ss(A3, B3, C3, 0);
Generator3d=c2d(Generator3, Ts);
A3d=Generator3d.A;
Generator4=ss(A4, B4, C4, 0);
Generator4d=c2d(Generator4, Ts);
A4d=Generator4d.A;
%  Overall system
AGS=zeros(nGS,nGS);
AGS(1:nTS,1:nTS)=A1d;
AGS(nTS+1:2*nTS,nTS+1:2*nTS)=A2d;
AGS(2*nTS+1:3*nTS,2*nTS+1:3*nTS)=A3d;
AGS(3*nTS+1:4*nTS,3*nTS+1:4*nTS)=A4d;
%  as in figure(1)
deltaP1=0;
deltaP2=0;
deltaP3=0;
deltaP4=0;
%   initial states
theta10=45*pi/180;    % 45 degrees
theta20=20*pi/180;
theta30=30*pi/180;
theta40=4*30*pi/180-(theta10+theta20+theta30);
x10=[theta10; 0; deltaP1/M1];
x20=[theta20; 0; deltaP2/M2];
x30=[theta30; 0; deltaP3/M3];
x40=[theta40; 0; deltaP4/M4];
X0=[x10; x20; x30; x40];   
%
X=[];                           % state as row vector
U=[];
YGS=[];
X(1,:)=X0;
Y=CGS*X0;                       % vector of current phases
YGS(1,:)=Y';                    % output as row vector
for k1=1:kend
    U=zeros(N1,1);
    for k2=1:N1
        for k3=1:N1
            U(k2,1)=U(k2,1)-AG(k2,k3)*sin(Y(k2)-Y(k3));
        end
    end
    X(k1+1,:)=AGS*X(k1,:)'+BGS*U;
    UGS(k1+1,:)=U';
    Y=CGS*X(k1+1,:)';
    YGS(k1+1,:)=Y';
end
subplot(2,1,1)
plot(Time,YGS(:,1:4),'b');
%axis([0 1 -0.5 0.5]);
% subplot(2,1,2)
% plot(Time,sin(YGS(:,1:2)+2*pi*50*Time'),'b')
% axis([0 0.35 -1.1 1.1]);

%%  Comparison with Kuramoto oscillator
%   for Mi=0
figure(5)
AK=0;
D1=80;  % in MW/Hz
D2=100;
D3=150;
D4=95;
BK1=1/D1;
BK2=1/D2;
BK3=1/D3;
BK4=1/D4;
CK=1;
DK=0;
GeneratorK1=ss(AK, BK1, CK, 0);
GeneratorK1d=c2d(GeneratorK1, Ts);
BK1d=GeneratorK1d.B;
GeneratorK2=ss(AK, BK2, CK, 0);
GeneratorK2d=c2d(GeneratorK2, Ts);
BK2d=GeneratorK2d.B;
GeneratorK3=ss(AK, BK3, CK, 0);
GeneratorK3d=c2d(GeneratorK3, Ts);
BK3d=GeneratorK3d.B;
GeneratorK4=ss(AK, BK4, CK, 0);
GeneratorK4d=c2d(GeneratorK4, Ts);
BK4d=GeneratorK4d.B;
AKbar=eye(4,4);
BKbar=zeros(4,4);
BKbar(1,1)=BK1d;
BKbar(2,2)=BK2d;
BKbar(3,3)=BK3d;
BKbar(4,4)=BK4d;
CKbar=eye(4,4);
theta10=45*pi/180;    % 45 degrees
theta20=20*pi/180;
theta30=30*pi/180;
theta40=4*30*pi/180-(theta10+theta20+theta30);
XK0=[theta10; theta20; theta30; theta40];
%
X=[];                           % state as row vector
UGS=[];
YGS=[];
X(1,:)=XK0';
Y=CKbar*XK0;                       % vector of current phases
YGS(1,:)=Y';                    % output as row vector
for k1=1:kend
    U=zeros(N1,1);
    for k2=1:N1
        for k3=1:N1
            U(k2,1)=U(k2,1)-AG(k2,k3)*sin(Y(k2)-Y(k3));
        end
    end
    X(k1+1,:)=AKbar*X(k1,:)'+BKbar*U;
    UGS(k1+1,:)=U';
    Y=CKbar*X(k1+1,:)';
    YGS(k1+1,:)=Y';
end
subplot(2,1,1)
plot(Time,YGS(:,1:4),'b');
%axis([0 1 -0.2 0.3]);
subplot(2,1,2)
plot(Time,sin(YGS(:,1:2)+2*pi*50*Time'),'b')
axis([0 0.15 -1.1 1.1]);



%%  use the linearised model:   sin(theta_i-theta_j) = theta_i-theta_j
%
figure(4)
Alin=zeros(nGS,nGS);
Alin(1:nTS,1:nTS)=A1;
Alin(nTS+1:2*nTS,nTS+1:2*nTS)=A2;
Alin(2*nTS+1:3*nTS,2*nTS+1:3*nTS)=A3;
Alin(3*nTS+1:4*nTS,3*nTS+1:4*nTS)=A4;
Blin=zeros(nGS,N1);
Blin(1:nTS,1)=B1;
Blin(nTS+1:2*nTS,2)=B2;
Blin(2*nTS+1:3*nTS,3)=B3;
Blin(3*nTS+1:4*nTS,4)=B4;
Clin=CGS;
L=diag(AG*ones(N1,1))-AG;
Alin=Alin-Blin*L*Clin;
Netlin=ss(Alin, Blin, Clin, 0);
Ylin=initial(Netlin, X0, Time);
plot(Time, Ylin);


%%  Disturbance by an additional load pL
%
figure(6)
pL=zeros(kend+1,4);
t1=2;     % begin of disturbance
t2=3;
pLMag=10; % 10 MW
for k1=ceil(t1/Ts):ceil(t2/Ts)
    pL(k1,1)=pLMag;
end
%
X=[];                           % state as row vector
U=[];
YGS=[];
X(1,:)=X0;
Y=CGS*X0;                       % vector of current phases
YGS(1,:)=Y';                    % output as row vector
for k1=1:kend
    U=zeros(N1,1);
    for k2=1:N1
        for k3=1:N1
            U(k2,1)=U(k2,1)-AG(k2,k3)*sin(Y(k2)-Y(k3));
        end
    end
    X(k1+1,:)=AGS*X(k1,:)'+BGS*U+BGS*pL(k1,:)';
    UGS(k1+1,:)=U';
    Y=CGS*X(k1+1,:)';
    YGS(k1+1,:)=Y';
end
subplot(2,1,1)
plot(Time,pL(:,1),'b');
subplot(2,1,2);
plot(Time,YGS(:,1:4),'b');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   Stability test for the four-area power system with individual dynamics
%
%  Agent modes with n=2
A1n=[0 1; 0 -D1/M1];
b1n=[0; 1/M1];
c1n=[1 0];
A2n=[0 1; 0 -D2/M2];
b2n=[0; 1/M2];
c2n=[1 0];
A3n=[0 1; 0 -D3/M3];
b3n=[0; 1/M3];
c3n=[1 0];
A4n=[0 1; 0 -D4/M4];
b4n=[0; 1/M4];
c4n=[1 0];
Ap1=-D1/M1;
Aq1hat=[1; 0];
bp1=1/M1;
bq1hat=[0; 0];
cp1=0;
%  graph
L=diag(sum(AG))-AG;
%  test matrix
Abreve=[Ap1 zeros(1,2) zeros(1,2) zeros(1,2);
        -Aq1hat A2n   zeros(2,2)  zeros(2,2);
        -Aq1hat zeros(2,2) A3n    zeros(2,2);
        -Aq1hat zeros(2,2)  zeros(2,2) A4n]...
       -[bp1 0 0 0;
         -bq1hat b2n zeros(2,1) zeros(2,1);
         -bq1hat zeros(2,1) b3n zeros(2,1);
         -bq1hat zeros(2,1) zeros(2,1) b4n]*L*...
         [cp1 zeros(1,2) zeros(1,2) zeros(1,2);
          0   c2n       zeros(1,2) zeros(1,2);
          0  zeros(1,2)   c3n  zeros(1,2);
          0  zeros(1,2)  zeros(1,2)  c4n];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Power system with partly renewable energy sources
%   Behaviour for balanced power generation
%
figure(7)
%  Areas 1 and 2
M1=2.5;   % inertia in s   instead of 5
M2=2;   %                instead of 5
D1=7.5*2*pi;  % in MW/Hz = 2pi MWs/rad
D2=7.2*2*pi;
%  new models
A1=[0    1   0;
    0 -D1/M1 1;
    0    0   0];
B1=[0; 1/M1; 0];
C1=[1 0 0];
A2=[0    1   0;
    0 -D2/M2 1;
    0    0   0];
B2=[0; 1/M2; 0];
C2=[1 0 0];
%   discrete-time models
Generator1=ss(A1, B1, C1, 0);
Generator1d=c2d(Generator1, Ts);
A1d=Generator1d.A;
B1d=Generator1d.B;
Generator2=ss(A2, B2, C2, 0);
Generator2d=c2d(Generator2, Ts);
A2d=Generator2d.A;
B2d=Generator2d.B;
%  Overall system
AGS=zeros(nGS,nGS);
AGS(1:nTS,1:nTS)=A1d;
AGS(nTS+1:2*nTS,nTS+1:2*nTS)=A2d;
AGS(2*nTS+1:3*nTS,2*nTS+1:3*nTS)=A3d;
AGS(3*nTS+1:4*nTS,3*nTS+1:4*nTS)=A4d;
BGS=zeros(nGS,N1);
BGS(1:nTS,1)=B1d;
BGS(nTS+1:2*nTS,2)=B2d;
BGS(2*nTS+1:3*nTS,3)=B3d;
BGS(3*nTS+1:4*nTS,4)=B4d;
CGS=zeros(N1, nGS);
CGS(1,1:nTS)=C1;
CGS(2,nTS+1:2*nTS)=C2;
CGS(3,2*nTS+1:3*nTS)=C3;
CGS(4,3*nTS+1:4*nTS)=C4;
%  initial states
deltaP1=0;
deltaP2=0;
deltaP3=0;
deltaP4=0;
theta1=0;
theta2=0;
theta3=0;
theta4=0;
%   initial states
theta10=45*pi/180;    % 45 degrees
theta20=20*pi/180;
theta30=30*pi/180;
theta40=4*30*pi/180-(theta10+theta20+theta30);
x10=[theta10; 0; deltaP1/M1];
x20=[theta20; 0; deltaP2/M2];
x30=[theta30; 0; deltaP3/M3];
x40=[theta40; 0; deltaP4/M4];
X0=[x10; x20; x30; x40];   
%
X=[];                           % state as row vector
U=[];
YGS=[];
X(1,:)=X0;
Y=CGS*X0;                       % vector of current phases
YGS(1,:)=Y';                    % output as row vector
for k1=1:kend
    U=zeros(N1,1);
    for k2=1:N1
        for k3=1:N1
            U(k2,1)=U(k2,1)-AG(k2,k3)*sin(Y(k2)-Y(k3));
        end
    end
    X(k1+1,:)=AGS*X(k1,:)'+BGS*U;
    UGS(k1+1,:)=U';
    Y=CGS*X(k1+1,:)';
    YGS(k1+1,:)=Y';
end
subplot(2,1,1)
plot(Time,YGS(:,1:4),'b');
axis([0 1 0.3 0.8]);
latexylabel('$$\tilde{\theta}_i$$ in rad $$\;$$');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 0.2 0.4 0.6 0.8 1]);
yticks([-0.2 0 0.2 0.4 0.6 0.8]);
latextitle('$${\tt SyncPowerSyst:}$$ Grid with renewable energy sources');
%
subplot(2,1,2)
plot(Time,110*sin(YGS(:,1:2)+2*pi*50*Time'),'b')
axis([0 0.15 -120 120]);
latexylabel('$$v_1, v_2$$ in kV');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 0.05 0.1 0.15]);
yticks([-110 0 110]);



%%   Disturbance by a load change
%
figure(8)
deltaP1=0;
deltaP2=0;
deltaP3=0;
deltaP4=0;
%  Disturbance by load change
pL=zeros(kend+1,4);
t1=0.2;     % begin of disturbance  0.2
t2=0.5;     % 0.5
pLMag=50; % 10 MW
for k1=ceil(t1/Ts):ceil(t2/Ts)
    pL(k1,1)=pLMag;
end
%   initial states
theta10=30*pi/180;   
theta20=30*pi/180; 
theta30=30*pi/180; 
theta40=30*pi/180; 
x10=[theta10; 0; deltaP1/M1];
x20=[theta20; 0; deltaP2/M2];
x30=[theta30; 0; deltaP3/M3];
x40=[theta40; 0; deltaP4/M4];
X0=[x10; x20; x30; x40];   
%
X=[];                           % state as row vector
U=[];
YGS=[];
X(1,:)=X0;
Y=CGS*X0;                       % vector of current phases
YGS(1,:)=Y';                    % output as row vector
for k1=1:kend
    U=zeros(N1,1);
    for k2=1:N1
        for k3=1:N1
            U(k2,1)=U(k2,1)-AG(k2,k3)*sin(Y(k2)-Y(k3));
        end
    end
    X(k1+1,:)=AGS*X(k1,:)'+BGS*U+BGS*pL(k1,:)';
    UGS(k1+1,:)=U';
    Y=CGS*X(k1+1,:)';
    YGS(k1+1,:)=Y';
end
subplot(2,1,1)
plot(Time,YGS(:,1:4),'b');
axis([0 1 0.4 0.7]);
latexylabel('$$\tilde{\theta}_i$$ in rad');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 0.2 0.4 0.6 0.8 1]);
yticks([-0.2 0 0.2 0.4 0.6 0.8]);
latextitle('$${\tt SyncPowerSyst:}$$ Grid with renewable energy sources');
%
subplot(2,1,2)
k1=ceil(0.1/Ts);
k2=ceil(0.2/Ts);
TimeV=Time(1,k1:k2)+0.65;   % shifted time
plot(Time, V1*sin(YGS(:,1:2)+2*pi*50*Time'),'b', TimeV, V1*sin(2*pi*50*TimeV'), 'b--');
axis([0.75 0.85 -120 120]);
latexylabel('$$v_1, v_2$$ in kV');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0.75 0.8 0.85]);
yticks([-110 0 110]);



%%  Figures
%
epsfigc15('SyncPowerSyst');

