Mirrors_Framework/Adaptive-PID-controller
1
0
Fork 0
mirror of https://github.com/NickNair/Adaptive-PID-controller.git synced 2026-02-05 09:49:42 +08:00
Code Issues Packages Projects Releases Wiki Activity

Added Matlab code

Browse Source
This commit is contained in:
Nikhil Nair
2022-06-08 22:59:59 +05:30
parent 0cfc443ab3
commit 2947cb2941
74 changed files with 2909 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

BIN
__pycache__/RBF.cpython-39.pyc
View File

Binary file not shown.

BIN
__pycache__/ac.cpython-39.pyc
View File

Binary file not shown.

BIN
__pycache__/neuralnetwork.cpython-39.pyc
View File

Binary file not shown.

BIN
__pycache__/singlearea.cpython-39.pyc
View File

Binary file not shown.

BIN
__pycache__/two_area.cpython-39.pyc
View File

Binary file not shown.

1
matlab_implementation/ACE1.txt Normal file
View File

File diff suppressed because one or more lines are too long

BIN
matlab_implementation/ACE1_comparison.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 27 KiB

1
matlab_implementation/ACE2.txt Normal file
View File

File diff suppressed because one or more lines are too long

BIN
matlab_implementation/ACE2_comparison.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 27 KiB

125
matlab_implementation/RBF.asv Normal file
View File

@@ -0,0 +1,125 @@
% Defining the class RBF
classdef RBF<handle
% Defining variables for the class
properties
X; % Variable for state
h; % Variable for number of hidden layers
mu; % Variable for mean
sigma; % Variable for variance
K; % Variable for K values of PID
Vprev; % Variable for previous value of value function
V; % Variable for value of value function
w; % Variable for Actor weight matrix
v; % Variable for Critic weight matrix
output; % Variable for hidden layer
aw; % Variable for learning parameter
av; % Variable for learning parameter
au; % Variable for learning parameter
asig; % Variable for learning parameter
gamma; % Variable for discount factor in RL
end
methods
% Defining the RBF constructor
function self = RBF(aw , av , au , asig , gamma,h )
self.X = zeros(3,1); % Initializing X
self.h = h; % Initializing for number of hidden layers h
self.mu = zeros(3,h); % Initializing mean (mu)
self.sigma = ones(1,h); % Initializing variance (sigma)
self.K = zeros(3,1); % Initializing K
self.Vprev = 0; % Initializing previous value function
self.V = 0; % Initializing value function
self.w = zeros(3,h); % Initializing w
self.v = zeros(1,h); % Initializing v
self.output = ones(h,1); % Initializing output
self.aw = aw; % Initializing learning parameter aw
self.av = av; % Initializing learning parameter av
self.au = au; % Initializing learning parameter au
self.asig = asig; % Initializing learning parameter asig
self.gamma = gamma; % Initializing learning parameter gamma
end
% Defining the function for the hidden layer of the RBF
function self = HiddenLayer(self)
% Looping over h hidden layers
for i = 1:self.h
% Gaussian Function (Equation 26 in the paper)
self.output(i) = exp(- (norm( self.X - self.mu(:,i) )^2 )/(2*self.sigma(i)^2) );
end
end
% Defining the function for output layer of the RBF
function self = OutputLayer(self)
% Equation 27 in the paper
self.K = self.w*self.output;
% Equating previous value function value to current value function value
self.Vprev = self.V;
% Equation 28 in the paper
self.V = self.v*self.output;
end
function self = Update(self, y_ref , Y, y1, y2 , y3 , delU )
% TD Error (Equation 30 in the paper)
del_TD = 0.5 * ( y_ref - Y )^2 + self.gamma*self.V - self.Vprev;
% Finding delY/delU
delY_delU = (Y - y1 )/delU ;
% Finding the gradient of the weight matrix
w_gradient = del_TD *delY_delU* ( [ (y1 - y2) ; -self.X(1) ; (y1 - 2*y2 + y3) ]*self.output');
% Update rule for w (Equation 32 in the paper)
self.w = self.w - self.aw * w_gradient;
% Finding the gradient of the mean (mu) matrix
mu_gradient = del_TD* (repmat(self.X,1,self.h) - self.mu).*( (( self.output.*(self.v') )./( (self.sigma).^2 )')' );
% Update rule for mean (Equation 34 in the paper)
self.mu = self.mu + self.au*mu_gradient;
temp = zeros(1,self.h);
for j = 1: self.h
temp(j) = norm( self.X - self.mu(:,j) )^2;
end
% Finding gradient for updating sigma (variance)
sig_gradient = (del_TD) * temp.* ((( self.output.*(self.v') )./( (self.sigma).^3 )')');
% Update rule for sigma (Equation 35 in the paper)
self.sigma = self.sigma + self.asig*sig_gradient;
% Update rule for critic weights, value function (Equation 33 in the paper)
self.v = self.v + self.av*del_TD*self.output';
end
% Defining function to initialize weights
function self = initialize_weights(self,area)
% Initialize the weights for area 1
if area == 1
self.w = readmatrix('w_matrix1.txt');
% Initialize the weights for area 2
elseif area == 2
self.w = readmatrix('w_matrix2.txt');
end
end
end
end

125
matlab_implementation/RBF.m Normal file
View File

@@ -0,0 +1,125 @@
% Defining the class RBF
classdef RBF<handle
% Defining variables for the class
properties
X; % Variable for state
h; % Variable for number of hidden layers
mu; % Variable for mean
sigma; % Variable for variance
K; % Variable for K values of PID
Vprev; % Variable for previous value of value function
V; % Variable for value of value function
w; % Variable for Actor weight matrix
v; % Variable for Critic weight matrix
output; % Variable for hidden layer
aw; % Variable for learning parameter
av; % Variable for learning parameter
au; % Variable for learning parameter
asig; % Variable for learning parameter
gamma; % Variable for discount factor in RL
end
methods
% Defining the RBF constructor
function self = RBF(aw , av , au , asig , gamma,h )
self.X = zeros(3,1); % Initializing X
self.h = h; % Initializing for number of hidden layers h
self.mu = zeros(3,h); % Initializing mean (mu)
self.sigma = ones(1,h); % Initializing variance (sigma)
self.K = zeros(3,1); % Initializing K
self.Vprev = 0; % Initializing previous value function
self.V = 0; % Initializing value function
self.w = zeros(3,h); % Initializing w
self.v = zeros(1,h); % Initializing v
self.output = ones(h,1); % Initializing output
self.aw = aw; % Initializing learning parameter aw
self.av = av; % Initializing learning parameter av
self.au = au; % Initializing learning parameter au
self.asig = asig; % Initializing learning parameter asig
self.gamma = gamma; % Initializing learning parameter gamma
end
% Defining the function for the hidden layer of the RBF
function self = HiddenLayer(self)
% Looping over h hidden layers
for i = 1:self.h
% Gaussian Function (Equation 26 in the paper)
self.output(i) = exp(- (norm( self.X - self.mu(:,i) )^2 )/(2*self.sigma(i)^2) );
end
end
% Defining the function for output layer of the RBF
function self = OutputLayer(self)
% Equation 27 in the paper
self.K = self.w*self.output;
% Equating previous value function value to current value function value
self.Vprev = self.V;
% Equation 28 in the paper
self.V = self.v*self.output;
end
function self = Update(self, y_ref , Y, y1, y2 , y3 , delU )
% TD Error (Equation 30 in the paper)
del_TD = 0.5 * ( y_ref - Y )^2 + self.gamma*self.V - self.Vprev;
% Finding delY/delU
delY_delU = (Y - y1 )/delU ;
% Finding the gradient of the weight matrix
w_gradient = del_TD *delY_delU* ( [ (y1 - y2) ; -self.X(1) ; (y1 - 2*y2 + y3) ]*self.output');
% Update rule for w (Equation 32 in the paper)
self.w = self.w - self.aw * w_gradient;
% Finding the gradient of the mean (mu) matrix
mu_gradient = del_TD* (repmat(self.X,1,self.h) - self.mu).*( (( self.output.*(self.v') )./( (self.sigma).^2 )')' );
% Update rule for mean (Equation 34 in the paper)
self.mu = self.mu + self.au*mu_gradient;
temp = zeros(1,self.h);
for j = 1: self.h
temp(j) = norm( self.X - self.mu(:,j) )^2;
end
% Finding gradient for updating sigma (variance)
sig_gradient = (del_TD) * temp.* ((( self.output.*(self.v') )./( (self.sigma).^3 )')');
% Update rule for sigma (Equation 35 in the paper)
self.sigma = self.sigma + self.asig*sig_gradient;
% Update rule for critic weights, value function (Equation 33 in the paper)
self.v = self.v + self.av*del_TD*self.output';
end
% Defining function to initialize weights
function self = initialize_weights(self,area)
% Initialize the weights for area 1
if area == 1
self.w = readmatrix('w_matrix1.txt');
% Initialize the weights for area 2
elseif area == 2
self.w = readmatrix('w_matrix2.txt');
end
end
end
end

116
matlab_implementation/TwoAreaPS.asv Normal file
View File

@@ -0,0 +1,116 @@
% Defining class TwoAreaPS
classdef TwoAreaPS<handle
% Defining the variables for the class
properties
yt_a1; % Defining variable yt for output of area 1
yt_a2; % Defining variable yt for output of area 2
Xprev; % Defining variable for previous value of X
X; % Defining variable X
Y; % Defining variable Y
Tg; % Defining variable for governor time constant
Tp; % Defining variable for plant time constant
Tt; % Defining variable for turbine time constant
Kp; % Defining variable for PID parameter Kp
T12; % Defining variable for tie line power constant
a12; % Defining variable for a12
R; % Defining variable for regulation
beta1; % Defining variable for frequency bias constant of area 1
beta2; % Defining variable for frequency bias constant of area 2
T; % Defining variable for time step
A; % Defining variable for state space model matrix A
B; % Defining variable for state space model matrix B
C; % Defining variable for state space model matrix C
Ad; % Defining variable for discretization parameter
Bd; % Defining variable for discretization parameter
end
methods
function self = TwoAreaPS( Tg, Tp, Tt, Kp, T12, a12, R, T, beta1, beta2, yt_a1, yt_a2)
self.yt_a1 = yt_a1; % Initializing variable for previous output value
self.yt_a2 = yt_a2; % Initializing variable for previous output value
self.Xprev = zeros( 7,1); % Initializing variable for input
self.X= zeros( 7,1); % Initializing variable for input
self.Y = zeros( 5,1 ); % Initializing variable for output
self.Tg = Tg; % Initializing variable for governor time constant
self.Tp = Tp ; % Initializing variable for plant time constant
self.Tt = Tt ; % Initializing variable for turbine time constant
self.Kp = Kp; % Initializing variable for PID parameter Kp
self.T12 = T12 ; % Initializing variable for tie line power constant
self.a12 = a12; % Initializing variable for a12
self.R = R ; % Initializing variable for regulation
self.beta1 = beta1; % Initializing variable for frequency bias constant of area 1
self.beta2 = beta2; % Initializing variable for frequency bias constant of area 2
self.T = T;
self.CalcDiscreteCoef();
end
function self = CalcDiscreteCoef(self)
Tg = self.Tg;
Tp = self.Tp;
Tt = self.Tt;
Kp = self.Kp;
T12 = self.T12 ;
a12 = self.a12;
R = self.R ;
% Defining state space model matrix A
self.A = [ -1/Tp Kp/Tp 0 -Kp/Tp 0 0 0 ;
0 -1/Tt 1/Tt 0 0 0 0 ;
-1/(R*Tg) 0 -1/Tg 0 0 0 0 ;
2*pi*T12 0 0 0 -2*pi*T12 0 0 ;
0 0 0 -Kp*a12/Tp -1/Tp Kp/Tp 0 ;
0 0 0 0 0 -1/Tt 1/Tt ;
0 0 0 0 -1/(R*Tg) 1/Tg -1/Tg ];
% Defining state space model matrix B
self.B = [ 0 0 -Kp/Tp 0 ;
0 0 0 0 ;
1/Tg 0 0 0 ;
0 0 0 0 ;
0 0 0 0 ;
0 0 0 0 ;
0 1/Tg 0 -Kp/Tp ] ;
% Defining state space model matrix C
self.C = [ self.beta1 0 0 1 0 0 0 ;
0 0 0 1 self.beta2 0 0 ;
0 0 0 1 0 0 0 ;
1 0 0 0 0 0 0 ;
0 0 0 0 1 0 0 ] ;
% Calculating system matrix for discrete time (eq 8 in the paper)
self.Ad = expm(self.A*self.T);
% Calculating input matrix for discrete time (eq 9 in the paper)
self.Bd = inv(self.A)*(self.Ad - eye(7) )*self.B;
end
% Defining function output
function self = Output(self,Ut)
% Updating previous values for area 1
self.yt_a1 = [self.Y(1) self.yt_a1(1) self.yt_a1(2)];
% Updating previous values for area 2
self.yt_a2 = [self.Y(2) self.yt_a2(1) self.yt_a2(2)];
% System equation (Equation 4 in the paper)
self.X = self.Ad*self.Xprev +self.Bd*Ut;
% Output equation (Equation 5 in the paper)
self.Y = self.C*self.Xprev;
% Updating previous values for X
self.Xprev = self.X;
end
end
end

116
matlab_implementation/TwoAreaPS.m Normal file
View File

@@ -0,0 +1,116 @@
% Defining class TwoAreaPS
classdef TwoAreaPS<handle
% Defining the variables for the class
properties
yt_a1; % Defining variable yt for output of area 1
yt_a2; % Defining variable yt for output of area 2
Xprev; % Defining variable for previous value of X
X; % Defining variable X
Y; % Defining variable Y
Tg; % Defining variable for governor time constant
Tp; % Defining variable for plant time constant
Tt; % Defining variable for turbine time constant
Kp; % Defining variable for PID parameter Kp
T12; % Defining variable for tie line power constant
a12; % Defining variable for a12
R; % Defining variable for regulation
beta1; % Defining variable for frequency bias constant of area 1
beta2; % Defining variable for frequency bias constant of area 2
T; % Defining variable for time step
A; % Defining variable for state space model matrix A
B; % Defining variable for state space model matrix B
C; % Defining variable for state space model matrix C
Ad; % Defining variable for discretization parameter
Bd; % Defining variable for discretization parameter
end
methods
function self = TwoAreaPS( Tg, Tp, Tt, Kp, T12, a12, R, T, beta1, beta2, yt_a1, yt_a2)
self.yt_a1 = yt_a1; % Initializing variable for previous output value
self.yt_a2 = yt_a2; % Initializing variable for previous output value
self.Xprev = zeros( 7,1); % Initializing variable for input
self.X= zeros( 7,1); % Initializing variable for input
self.Y = zeros( 5,1 ); % Initializing variable for output
self.Tg = Tg; % Initializing variable for governor time constant
self.Tp = Tp ; % Initializing variable for plant time constant
self.Tt = Tt ; % Initializing variable for turbine time constant
self.Kp = Kp; % Initializing variable for PID parameter Kp
self.T12 = T12 ; % Initializing variable for tie line power constant
self.a12 = a12; % Initializing variable for a12
self.R = R ; % Initializing variable for regulation
self.beta1 = beta1; % Initializing variable for frequency bias constant of area 1
self.beta2 = beta2; % Initializing variable for frequency bias constant of area 2
self.T = T;
self.CalcDiscreteCoef();
end
function self = CalcDiscreteCoef(self)
Tg = self.Tg;
Tp = self.Tp;
Tt = self.Tt;
Kp = self.Kp;
T12 = self.T12 ;
a12 = self.a12;
R = self.R ;
% Defining state space model matrix A
self.A = [ -1/Tp Kp/Tp 0 -Kp/Tp 0 0 0 ;
0 -1/Tt 1/Tt 0 0 0 0 ;
-1/(R*Tg) 0 -1/Tg 0 0 0 0 ;
2*pi*T12 0 0 0 -2*pi*T12 0 0 ;
0 0 0 -Kp*a12/Tp -1/Tp Kp/Tp 0 ;
0 0 0 0 0 -1/Tt 1/Tt ;
0 0 0 0 -1/(R*Tg) 1/Tg -1/Tg ];
% Defining state space model matrix B
self.B = [ 0 0 -Kp/Tp 0 ;
0 0 0 0 ;
1/Tg 0 0 0 ;
0 0 0 0 ;
0 0 0 0 ;
0 0 0 0 ;
0 1/Tg 0 -Kp/Tp ] ;
% Defining state space model matrix C
self.C = [ self.beta1 0 0 1 0 0 0 ;
0 0 0 1 self.beta2 0 0 ;
0 0 0 1 0 0 0 ;
1 0 0 0 0 0 0 ;
0 0 0 0 1 0 0 ] ;
% Calculating system matrix for discrete time (eq 8 in the paper)
self.Ad = expm(self.A*self.T);
% Calculating input matrix for discrete time (eq 9 in the paper)
self.Bd = inv(self.A)*(self.Ad - eye(7) )*self.B;
end
% Defining function output
function self = Output(self,Ut)
% Updating previous values for area 1
self.yt_a1 = [self.Y(1) self.yt_a1(1) self.yt_a1(2)];
% Updating previous values for area 2
self.yt_a2 = [self.Y(2) self.yt_a2(1) self.yt_a2(2)];
% System equation (Equation 4 in the paper)
self.X = self.Ad*self.Xprev +self.Bd*Ut;
% Output equation (Equation 5 in the paper)
self.Y = self.C*self.Xprev;
% Updating previous values for X
self.Xprev = self.X;
end
end
end

162
matlab_implementation/TwoAreaPS_PID.asv Normal file
View File

@@ -0,0 +1,162 @@
% Load Disturbance for area 1
PL1 = 0;
% Load Disturbance for area 2
PL2 = 0;
% Defining and initializing variables
ut1 = []; % Variable for control signal for area 1
ut2 = []; % Variable for control signal for area 2
pl1 = []; % Variable for load disturbance for area 1
pl2 = []; % Variable for load disturbance for area 2
ACE1 = []; % Variable for area control error for area 1
ACE2 = []; % Variable for area control error for area 2
TieLine = []; % Variable for TieLine Power
Freq1 =[]; % Variable for Frequency change for area 2
Freq2 = [];% Variable for Frequency change for area 2
time = []; % Variable for time
Tg = 0.08; % Initializing variable for governor time constant
Tp = 20; % Initializing variable for plant time constant
Tt = 0.3; % Initializing variable for Turbine time constant
Kp = 120; % Initliazing PID parameter K
T12 = 0.545/(2*pi); % Initializing variable for Tie line power constant
a12 = -1;
R = 5; % Initializing variable for Regulation
T = 1/80;
dt = 1/80;
beta1 = 0.425; % Initializing variable for Frequency bias constant of area 1
beta2 = 0.425; % Initializing variable for Frequency bias constant of area 2
% Defining object System for the class TwoAreaPS
System = TwoAreaPS( Tg, Tp, Tt, Kp, T12, a12, R, T, beta1, beta2, [0,0,0],[0,0,0] );
KP1 = 1.723845;
KI1 = 0.014555;
KD1 = 17.936909;
KP2 = 0.106744;
KI2 = 0.006323;
KD2 = 15.023620;
ut_1 = 0;
ut_2 = 0;
t = 200;
y = [];
x = [];
for i = 1:t/dt
% Finding error signal e_t for area 1
e_t_a1 = 0 - System.yt_a1(1);
% Finding derivative of the output (delta y) for area 1
del_y_a1 = System.yt_a1(1) - System.yt_a1(2);
% Finding double derivate of the output (delta square y) for area 1
del2_y_a1 = System.yt_a1(1) - 2*System.yt_a1(2) + System.yt_a1(3);
% Controller definition (Equation 22 in the paper) for area 1
ut_1 = ut_1 + (KI1*e_t_a1 - KP1*del_y_a1 - KD1*del2_y_a1);
% Controller definition (Equation 22 in the paper) for area 2
ut_2 = ut_2 + (KI2*e_t_a2 - KP2*del_y_a2 - KD2*del2_y_a2);
% Defining control signal ut1 for area 1
ut1 = [ut1 ut_1];
% Defining control signal ut2 for area 2
ut2 = [ut2 ut_2];
if (i*dt >=40 && i*dt <= 80)
PL1 = 0.003;
elseif( i*dt > 80 && i*dt <= 120)
PL1 = 0.006;
elseif( i*dt > 120 && i*dt <= 160)
PL1 = 0.009;
else
PL1 = 0;
end
PL2 = 0;
pl1 = [pl1 PL1];
pl2 = [pl2 PL2];
% Defining control signal Ut
Ut = [ ut_1 ; ut_2 ; PL1 ; PL2 ];
System.Output(Ut);
% Defining ACE for area 1
ACE1 = [ACE1 System.Y(1)];
% Defining ACE for area 2
ACE2 = [ACE2 System.Y(2)];
% Defining Tie Line Power
TieLine = [ TieLine System.Y(3)];
% Defining Frequency Change for area 1
Freq1 = [ Freq1 System.Y(4)];
% Defining Frequency Change for area 2
Freq2 = [ Freq2 System.Y(5)];
% Defining variable time for our time period
time = [time i*dt];
ITAE = [];
ITAE = [ITAE, dt*(i*dt)*( abs(System.Y(1)) + abs(System.Y(2)))];
ISTE = [];
ISTE = [ISTE, dt*(i*dt)*( System.Y(1)^2 + System.Y(2)^2 )];
IAE = [];
IAE = [IAE, dt*( abs(System.Y(1)) + abs(System.Y(2)) ) ];
ISE = [];
ISE = [ISE, dt*( System.Y(1)^2 + System.Y(2)^2 ) ];
end
fprintf('IAE is');
disp(sum(IAE));
fprintf('ISE is');
disp(sum(ISE));
fprintf('ISTE is');
disp(sum(ISTE));
fprintf('ITAE is');
disp(sum(ITAE));
figure(1);
subplot(2,1,1);
plot(time,pl1 );
title(' Area 1 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
subplot(2,1,2);
plot(time, pl2);
title('Area 2 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
figure(2);
subplot(2,1,1);
plot(time,ut1 );
title(' Area 1 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
subplot(2,1,2);
plot(time, ut2);
title('Area 2 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
figure(3);
subplot(2,1,1);
plot(time, ACE1);
title(' Area 1 ACE vs Time for PSO Tuned PID')
xlabel('Time(s)')
ylabel('ACE (pu)')
subplot(2,1,2);
plot(time, ACE2);
title('Area 2 ACE vs Time for PSO Tuned PID')
xlabel('Time(s)')
ylabel('ACE (pu)')
ACE1_pso = ACE1;
ACE2_pso = ACE2;

162
matlab_implementation/TwoAreaPS_PID.m Normal file
View File

@@ -0,0 +1,162 @@
% Load Disturbance for area 1
PL1 = 0;
% Load Disturbance for area 2
PL2 = 0;
% Defining and initializing variables
ut1 = []; % Variable for control signal for area 1
ut2 = []; % Variable for control signal for area 2
pl1 = []; % Variable for load disturbance for area 1
pl2 = []; % Variable for load disturbance for area 2
ACE1 = []; % Variable for area control error for area 1
ACE2 = []; % Variable for area control error for area 2
TieLine = []; % Variable for TieLine Power
Freq1 =[]; % Variable for Frequency change for area 2
Freq2 = [];% Variable for Frequency change for area 2
time = []; % Variable for time
Tg = 0.08; % Initializing variable for governor time constant
Tp = 20; % Initializing variable for plant time constant
Tt = 0.3; % Initializing variable for Turbine time constant
Kp = 120; % Initliazing PID parameter K
T12 = 0.545/(2*pi); % Initializing variable for Tie line power constant
a12 = -1;
R = 5; % Initializing variable for Regulation
T = 1/80;
dt = 1/80;
beta1 = 0.425; % Initializing variable for Frequency bias constant of area 1
beta2 = 0.425; % Initializing variable for Frequency bias constant of area 2
% Defining object System for the class TwoAreaPS
System = TwoAreaPS( Tg, Tp, Tt, Kp, T12, a12, R, T, beta1, beta2, [0,0,0],[0,0,0] );
KP1 = 1.723845;
KI1 = 0.014555;
KD1 = 17.936909;
KP2 = 0.106744;
KI2 = 0.006323;
KD2 = 15.023620;
ut_1 = 0;
ut_2 = 0;
t = 200;
y = [];
x = [];
for i = 1:t/dt
% Finding error signal e_t for area 1
e_t_a1 = 0 - System.yt_a1(1);
% Finding derivative of the output (delta y) for area 1
del_y_a1 = System.yt_a1(1) - System.yt_a1(2);
% Finding double derivate of the output (delta square y) for area 1
del2_y_a1 = System.yt_a1(1) - 2*System.yt_a1(2) + System.yt_a1(3);
% Controller definition (Equation 22 in the paper) for area 1
ut_1 = ut_1 + (KI1*e_t_a1 - KP1*del_y_a1 - KD1*del2_y_a1);
% Controller definition (Equation 22 in the paper) for area 2
ut_2 = ut_2 + (KI2*e_t_a2 - KP2*del_y_a2 - KD2*del2_y_a2);
% Defining control signal ut1 for area 1
ut1 = [ut1 ut_1];
% Defining control signal ut2 for area 2
ut2 = [ut2 ut_2];
if (i*dt >=40 && i*dt <= 80)
PL1 = 0.003;
elseif( i*dt > 80 && i*dt <= 120)
PL1 = 0.006;
elseif( i*dt > 120 && i*dt <= 160)
PL1 = 0.009;
else
PL1 = 0;
end
PL2 = 0;
pl1 = [pl1 PL1];
pl2 = [pl2 PL2];
% Defining control signal Ut
Ut = [ ut_1 ; ut_2 ; PL1 ; PL2 ];
System.Output(Ut);
% Defining ACE for area 1
ACE1 = [ACE1 System.Y(1)];
% Defining ACE for area 2
ACE2 = [ACE2 System.Y(2)];
% Defining Tie Line Power
TieLine = [ TieLine System.Y(3)];
% Defining Frequency Change for area 1
Freq1 = [ Freq1 System.Y(4)];
% Defining Frequency Change for area 2
Freq2 = [ Freq2 System.Y(5)];
% Defining variable time for our time period
time = [time i*dt];
ITAE = [];
ITAE = [ITAE, dt*(i*dt)*( abs(System.Y(1)) + abs(System.Y(2)))];
ISTE = [];
ISTE = [ISTE, dt*(i*dt)*( System.Y(1)^2 + System.Y(2)^2 )];
IAE = [];
IAE = [IAE, dt*( abs(System.Y(1)) + abs(System.Y(2)) ) ];
ISE = [];
ISE = [ISE, dt*( System.Y(1)^2 + System.Y(2)^2 ) ];
end
fprintf('IAE is');
disp(sum(IAE));
fprintf('ISE is');
disp(sum(ISE));
fprintf('ISTE is');
disp(sum(ISTE));
fprintf('ITAE is');
disp(sum(ITAE));
figure(1);
subplot(2,1,1);
plot(time,pl1 );
title(' Area 1 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
subplot(2,1,2);
plot(time, pl2);
title('Area 2 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
figure(2);
subplot(2,1,1);
plot(time,ut1 );
title(' Area 1 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
subplot(2,1,2);
plot(time, ut2);
title('Area 2 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
figure(3);
subplot(2,1,1);
plot(time, ACE1);
title(' Area 1 ACE vs Time for PSO Tuned PID')
xlabel('Time(s)')
ylabel('ACE (pu)')
subplot(2,1,2);
plot(time, ACE2);
title('Area 2 ACE vs Time for PSO Tuned PID')
xlabel('Time(s)')
ylabel('ACE (pu)')
ACE1_pso = ACE1;
ACE2_pso = ACE2;

192
matlab_implementation/TwoAreaPS_PID_2.asv Normal file
View File

@@ -0,0 +1,192 @@
% Load Disturbance for area 1
PL1 = 0;
% Load Disturbance for area 2
PL2 = 0;
% Defining and initializing variables
ut1 = []; % Variable for control signal for area 1
ut2 = []; % Variable for control signal for area 2
pl1 = []; % Variable for load disturbance for area 1
pl2 = []; % Variable for load disturbance for area 2
ACE1 = []; % Variable for area control error for area 1
ACE2 = []; % Variable for area control error for area 2
TieLine = []; % Variable for TieLine Power
Freq1 =[]; % Variable for Frequency change for area 2
Freq2 = [];% Variable for Frequency change for area 2
time = []; % Variable for time
Tg = 0.08; % Initializing variable for governor time constant
Tp = 20; % Initializing variable for plant time constant
Tt = 0.3; % Initializing variable for Turbine time constant
Kp = 120; % Initliazing PID parameter K
T12 = 0.545/(2*pi); % Initializing variable for Tie line power constant
a12 = -1;
R = 5; % Initializing variable for Regulation
T = 1/80;
dt = 1/80;
beta1 = 0.425; % Initializing variable for Frequency bias constant of area 1
beta2 = 0.425; % Initializing variable for Frequency bias constant of area 2
% Defining object System for the class TwoAreaPS
System = TwoAreaPS( Tg, Tp, Tt, Kp, T12, a12, R, T, beta1, beta2, [0,0,0],[0,0,0] );
KP1 = 1.723845;
KI1 = 0.014555;
KD1 = 17.936909;
KP2 = 0.106744;
KI2 = 0.006323;
KD2 = 15.023620;
ut_1 = 0;
ut_2 = 0;
t = 200;
y = [];
x = [];
for i = 1:t/dt
% Finding error signal e_t for area 1
e_t_a1 = 0 - System.yt_a1(1);
% Finding derivative of the output (delta y) for area 1
del_y_a1 = System.yt_a1(1) - System.yt_a1(2);
% Finding double derivate of the output (delta square y) for area 1
del2_y_a1 = System.yt_a1(1) - 2*System.yt_a1(2) + System.yt_a1(3);
e_t_a2 = 0 - System.yt_a2(1);
% Finding derivative of the output (delta y) for area 1
del_y_a2 = System.yt_a2(1) - System.yt_a2(2);
% Finding double derivate of the output (delta square y) for area 1
del2_y_a2 = System.yt_a2(1) - 2*System.yt_a2(2) + System.yt_a2(3);
% Controller definition (Equation 22 in the paper) for area 1
ut_1 = ut_1 + (KI1*e_t_a1 - KP1*del_y_a1 - KD1*del2_y_a1);
% Controller definition (Equation 22 in the paper) for area 2
ut_2 = ut_2 + (KI2*e_t_a2 - KP2*del_y_a2 - KD2*del2_y_a2);
% Defining control signal ut1 for area 1
ut1 = [ut1 ut_1];
% Defining control signal ut2 for area 2
ut2 = [ut2 ut_2];
if (i*dt >=40 && i*dt <= 80)
PL1 = 0.003;
elseif( i*dt > 80 && i*dt <= 120)
PL1 = 0.006;
elseif( i*dt > 120 && i*dt <= 160)
PL1 = 0.009;
else
PL1 = 0;
end
PL2 = 0;
pl1 = [pl1 PL1];
pl2 = [pl2 PL2];
% Defining control signal Ut
Ut = [ ut_1 ; ut_2 ; PL1 ; PL2 ];
System.Output(Ut);
% Defining ACE for area 1
ACE1 = [ACE1 System.Y(1)];
% Defining ACE for area 2
ACE2 = [ACE2 System.Y(2)];
% Defining Tie Line Power
TieLine = [ TieLine System.Y(3)];
% Defining Frequency Change for area 1
Freq1 = [ Freq1 System.Y(4)];
% Defining Frequency Change for area 2
Freq2 = [ Freq2 System.Y(5)];
% Defining variable time for our time period
time = [time i*dt];
ITAE = [];
ITAE = [ITAE, dt*(i*dt)*( abs(System.Y(1)) + abs(System.Y(2)))];
ISTE = [];
ISTE = [ISTE, dt*(i*dt)*( System.Y(1)^2 + System.Y(2)^2 )];
IAE = [];
IAE = [IAE, dt*( abs(System.Y(1)) + abs(System.Y(2)) ) ];
ISE = [];
ISE = [ISE, dt*( System.Y(1)^2 + System.Y(2)^2 ) ];
end
fprintf('IAE is');
disp(sum(IAE));
fprintf('ISE is');
disp(sum(ISE));
fprintf('ISTE is');
disp(sum(ISTE));
fprintf('ITAE is');
disp(sum(ITAE));
figure(1);
subplot(2,1,1);
plot(time,pl1 );
title(' Area 1 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
subplot(2,1,2);
plot(time, pl2);
title('Area 2 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
figure(2);
subplot(2,1,1);
plot(time,ut1 );
title(' Area 1 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
subplot(2,1,2);
plot(time, ut2);
title('Area 2 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
figure(3);
subplot(2,1,1);
plot(time, ACE1);
title(' Area 1 ACE vs Time')
xlabel('Time(s)')
ylabel('ACE (pu)')
subplot(2,1,2);
plot(time, ACE2);
title('Area 2 ACE vs Time')
xlabel('Time(s)')
ylabel('ACE (pu)')
figure(4);
subplot(2,1,1);
plot(time, Freq1);
title(' Area 1 Frequency Change vs Time')
xlabel('Time(s)')
ylabel('Frequency Change (pu)')
subplot(2,1,2);
plot(time, Freq2);
title('Area 2 Frequency Change vs Time')
xlabel('Time(s)')
ylabel('Frequency Change (pu)')
figure(5);
plot(time, TieLine);
title('Tie Line Power vs Time')
xlabel('Time(s)')
ylabel('Tie Line Power (pu)')
ACE1_pso = ACE1;
ACE2_pso = ACE2;
pso_freq1 = Freq1;
pso_freq2 = Freq2;
pso_tileine = TieLine;
pso_ut1= ut1 ;
pso_ut2 = ut2 ;

192
matlab_implementation/TwoAreaPS_PID_2.m Normal file
View File

@@ -0,0 +1,192 @@
% Load Disturbance for area 1
PL1 = 0;
% Load Disturbance for area 2
PL2 = 0;
% Defining and initializing variables
ut1 = []; % Variable for control signal for area 1
ut2 = []; % Variable for control signal for area 2
pl1 = []; % Variable for load disturbance for area 1
pl2 = []; % Variable for load disturbance for area 2
ACE1 = []; % Variable for area control error for area 1
ACE2 = []; % Variable for area control error for area 2
TieLine = []; % Variable for TieLine Power
Freq1 =[]; % Variable for Frequency change for area 2
Freq2 = [];% Variable for Frequency change for area 2
time = []; % Variable for time
Tg = 0.08; % Initializing variable for governor time constant
Tp = 20; % Initializing variable for plant time constant
Tt = 0.3; % Initializing variable for Turbine time constant
Kp = 120; % Initliazing PID parameter K
T12 = 0.545/(2*pi); % Initializing variable for Tie line power constant
a12 = -1;
R = 5; % Initializing variable for Regulation
T = 1/80;
dt = 1/80;
beta1 = 0.425; % Initializing variable for Frequency bias constant of area 1
beta2 = 0.425; % Initializing variable for Frequency bias constant of area 2
% Defining object System for the class TwoAreaPS
System = TwoAreaPS( Tg, Tp, Tt, Kp, T12, a12, R, T, beta1, beta2, [0,0,0],[0,0,0] );
KP1 = 1.723845;
KI1 = 0.014555;
KD1 = 17.936909;
KP2 = 0.106744;
KI2 = 0.006323;
KD2 = 15.023620;
ut_1 = 0;
ut_2 = 0;
t = 200;
y = [];
x = [];
for i = 1:t/dt
% Finding error signal e_t for area 1
e_t_a1 = 0 - System.yt_a1(1);
% Finding derivative of the output (delta y) for area 1
del_y_a1 = System.yt_a1(1) - System.yt_a1(2);
% Finding double derivate of the output (delta square y) for area 1
del2_y_a1 = System.yt_a1(1) - 2*System.yt_a1(2) + System.yt_a1(3);
e_t_a2 = 0 - System.yt_a2(1);
% Finding derivative of the output (delta y) for area 1
del_y_a2 = System.yt_a2(1) - System.yt_a2(2);
% Finding double derivate of the output (delta square y) for area 1
del2_y_a2 = System.yt_a2(1) - 2*System.yt_a2(2) + System.yt_a2(3);
% Controller definition (Equation 22 in the paper) for area 1
ut_1 = ut_1 + (KI1*e_t_a1 - KP1*del_y_a1 - KD1*del2_y_a1);
% Controller definition (Equation 22 in the paper) for area 2
ut_2 = ut_2 + (KI2*e_t_a2 - KP2*del_y_a2 - KD2*del2_y_a2);
% Defining control signal ut1 for area 1
ut1 = [ut1 ut_1];
% Defining control signal ut2 for area 2
ut2 = [ut2 ut_2];
if (i*dt >=40 && i*dt <= 80)
PL1 = 0.003;
elseif( i*dt > 80 && i*dt <= 120)
PL1 = 0.006;
elseif( i*dt > 120 && i*dt <= 160)
PL1 = 0.009;
else
PL1 = 0;
end
PL2 = 0;
pl1 = [pl1 PL1];
pl2 = [pl2 PL2];
% Defining control signal Ut
Ut = [ ut_1 ; ut_2 ; PL1 ; PL2 ];
System.Output(Ut);
% Defining ACE for area 1
ACE1 = [ACE1 System.Y(1)];
% Defining ACE for area 2
ACE2 = [ACE2 System.Y(2)];
% Defining Tie Line Power
TieLine = [ TieLine System.Y(3)];
% Defining Frequency Change for area 1
Freq1 = [ Freq1 System.Y(4)];
% Defining Frequency Change for area 2
Freq2 = [ Freq2 System.Y(5)];
% Defining variable time for our time period
time = [time i*dt];
ITAE = [];
ITAE = [ITAE, dt*(i*dt)*( abs(System.Y(1)) + abs(System.Y(2)))];
ISTE = [];
ISTE = [ISTE, dt*(i*dt)*( System.Y(1)^2 + System.Y(2)^2 )];
IAE = [];
IAE = [IAE, dt*( abs(System.Y(1)) + abs(System.Y(2)) ) ];
ISE = [];
ISE = [ISE, dt*( System.Y(1)^2 + System.Y(2)^2 ) ];
end
fprintf('IAE is');
disp(sum(IAE));
fprintf('ISE is');
disp(sum(ISE));
fprintf('ISTE is');
disp(sum(ISTE));
fprintf('ITAE is');
disp(sum(ITAE));
figure(1);
subplot(2,1,1);
plot(time,pl1 );
title(' Area 1 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
subplot(2,1,2);
plot(time, pl2);
title('Area 2 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
figure(2);
subplot(2,1,1);
plot(time,ut1 );
title(' Area 1 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
subplot(2,1,2);
plot(time, ut2);
title('Area 2 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
figure(3);
subplot(2,1,1);
plot(time, ACE1);
title(' Area 1 ACE vs Time')
xlabel('Time(s)')
ylabel('ACE (pu)')
subplot(2,1,2);
plot(time, ACE2);
title('Area 2 ACE vs Time')
xlabel('Time(s)')
ylabel('ACE (pu)')
figure(4);
subplot(2,1,1);
plot(time, Freq1);
title(' Area 1 Frequency Change vs Time')
xlabel('Time(s)')
ylabel('Frequency Change (pu)')
subplot(2,1,2);
plot(time, Freq2);
title('Area 2 Frequency Change vs Time')
xlabel('Time(s)')
ylabel('Frequency Change (pu)')
figure(5);
plot(time, TieLine);
title('Tie Line Power vs Time')
xlabel('Time(s)')
ylabel('Tie Line Power (pu)')
ACE1_pso = ACE1;
ACE2_pso = ACE2;
pso_freq1 = Freq1;
pso_freq2 = Freq2;
pso_tileine = TieLine;
pso_ut1= ut1 ;
pso_ut2 = ut2 ;

224
matlab_implementation/TwoArea_Adaptive.m Normal file
View File

@@ -0,0 +1,224 @@
% Load disturbance for Area 1
PL1 = 0;
% Load disturbance for Area 2
PL2 = 0;
% Importing RBF Class
import RBF.*
% Defining an object (nn1) of the class RBF
nn1 = RBF( 8e+1, 33.3e-2, 5e+1 , 9e-1, 0.95, 3);
% Defining an object (nn2) of the class RBF
nn2 = RBF( 8e+1, 33.3e-2, 5e+1 , 9e-1, 0.95, 3);
% Defining and initializing variables
ut1 = []; % Variable for control signal for area 1
ut2 = []; % Variable for control signal for area 2
pl1 = []; % Variable for load disturbance for area 1
pl2 = []; % Variable for load disturbance for area 2
ACE1 = []; % Variable for area control error for area 1
ACE2 = []; % Variable for area control error for area 2
TieLine = []; % Variable for TieLine Power
Freq1 =[]; % Variable for Frequency change for area 2
Freq2 = [];% Variable for Frequency change for area 2
time = []; % Variable for time
Kp_data = [];
Tg = 0.08; % Initializing variable for governor time constant
Tp = 20; % Initializing variable for plant time constant
Tt = 0.3; % Initializing variable for Turbine time constant
Kp = 120; % Initializing variable for PID parameter Kp
T12 = 0.545/(2*pi); % Initializing variable for Tie line power constant
a12 = -1; % Initializing variable a12
R = 5; % Initializing variable for Regulation
T = 1/80; % Initializing Time step
dt = 1/80; % Initializing Time step
beta1 = 0.425; % Initializing variable for Frequency bias constant of area 1
beta2 = 0.425; % Initializing variable for Frequency bias constant of area 2
% Defining object System for the class TwoAreaPS
System = TwoAreaPS( Tg, Tp, Tt, Kp, T12, a12, R, T, beta1, beta2, [0,0,0],[0,0,0] );
% Initializing Ki
Ki = 0 ;
% Initializing Kd
Kd = 0 ;
% Initializing Kp
Kp = 0 ;
% Initializing weights
nn1.initialize_weights(1);
nn2.initialize_weights(2);
% Initializing control signals
ut_1 = 0;
ut_2 = 0;
% Duration of run
t = 200;
% Output variable
y = [];
for i = 1:t/dt
% Finding error signal e_t for area 1
e_t_a1 = 0 - System.yt_a1(1);
% Finding derivative of the output (delta y) for area 1
del_y_a1 = System.yt_a1(1) - System.yt_a1(2);
% Finding double derivate of the output (delta square y) for area 1
del2_y_a1 = System.yt_a1(1) - 2*System.yt_a1(2) + System.yt_a1(3);
% Defining the state space X for area 1
nn1.X = [ e_t_a1 ; -del_y_a1 ; -del2_y_a1];
% Calling HiddenLayer function using object nn1
nn1.HiddenLayer();
% Calling OutputLayer function using object nn1
nn1.OutputLayer();
% Finding error signal e_t for area 2
e_t_a2 = 0 - System.yt_a2(1);
% Finding derivative of the output (delta y) for area 2
del_y_a2 = System.yt_a2(1) - System.yt_a2(2);
% Finding double derivate of the output (delta square y) for area 2
del2_y_a2 = System.yt_a2(1) - 2*System.yt_a2(2) + System.yt_a2(3);
% % Defining the state space X for area 2
nn2.X = [ e_t_a2 ; -del_y_a2 ; -del2_y_a2];
% Calling HiddenLayer function using object nn2
nn2.HiddenLayer();
% Calling OutputLayer function using object nn2
nn2.OutputLayer();
% Controller definition (Equation 22 in the paper) for area 1
ut_1 = ut_1 + (nn1.K(2)*e_t_a1 - nn1.K(1)*del_y_a1 - nn1.K(3)*del2_y_a1);
% Controller definition (Equation 22 in the paper) for area 2
ut_2 = ut_2 + (nn1.K(2)*e_t_a2 - nn1.K(1)*del_y_a2 - nn1.K(3)*del2_y_a2);
% Defining control signal ut1 for area 1
ut1 = [ut1 ut_1];
% Defining control signal ut2 for area 2
ut2 = [ut2 ut_2];
% Load disturbances
if (i*dt >=40 && i*dt <= 80)
PL1 = 0.003;
elseif( i*dt > 80 && i*dt <= 120)
PL1 = 0.006;
elseif( i*dt > 120 && i*dt <= 160)
PL1 = 0.009;
else
PL1 = 0;
end
PL2 = 0;
pl1 = [pl1 PL1];
pl2 = [pl2 PL2];
% Defining control signal Ut
Ut = [ ut_1 ; ut_2 ; PL1 ; PL2 ];
System.Output(Ut);
% Calling Update function using the object nn1 for area 1
nn1.Update(0 ,System.Y(1), System.yt_a1(1), System.yt_a1(2), System.yt_a1(3));
% Calling Update function using the object nn2 for area 2
nn2.Update(0 ,System.Y(2), System.yt_a2(1), System.yt_a2(2), System.yt_a2(3) );
% Defining ACE for area 1
ACE1 = [ACE1 System.Y(1)];
% Defining ACE for area 2
ACE2 = [ACE2 System.Y(2)];
% Defining Tie Line Power
TieLine = [ TieLine System.Y(3)];
% Defining Frequency Change for area 1
Freq1 = [ Freq1 System.Y(4)];
% Defining Frequency Change for area 2
Freq2 = [ Freq2 System.Y(5)];
% Defining variable time for our time period
time = [time i*dt];
% Metrics calculations
%Integral time absolute error(eq 38 in paper)
ITAE = [];
ITAE = [ITAE, dt*(i*dt)*( abs(System.Y(1)) + abs(System.Y(2)))];
% Integral time squared error(eq 36 in paper)
ITSE = [];
ITSE = [ITSE, dt*(i*dt)*( System.Y(1)^2 + System.Y(2)^2 )];
% Integral absolute error(eq 37 in paper)
IAE = [];
IAE = [IAE, dt*( abs(System.Y(1)) + abs(System.Y(2)) ) ];
% Integral Square error(eq 39 in paper)
ISE = [];
ISE = [ISE, dt*( System.Y(1)^2 + System.Y(2)^2 ) ];
end
fprintf('IAE is');
disp(sum(IAE));
fprintf('ISE is');
disp(sum(ISE));
fprintf('ITSE is');
disp(sum(ITSE));
fprintf('ITAE is');
disp(sum(ITAE));
figure(1);
subplot(2,1,1);
plot(time,pl1 );
title(' Area 1 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
subplot(2,1,2);
plot(time, pl2);
title('Area 2 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
figure(2);
subplot(2,1,1);
plot(time,ut1 );
title(' Area 1 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
subplot(2,1,2);
plot(time, ut2);
title('Area 2 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
% Plotting ACE vs Time
figure(3);
subplot(2,1,1);
plot(time, ACE1);
title(' Area 1 ACE vs Time')
xlabel('Time(s)')
ylabel('ACE (pu)')
subplot(2,1,2);
plot(time, ACE2);
title('Area 2 ACE vs Time')
xlabel('Time(s)')
ylabel('ACE (pu)')
figure(4);
subplot(2,1,1);
plot(time, Freq1);
title(' Area 1 Frequency Change vs Time')
xlabel('Time(s)')
ylabel('Frequency Change (pu)')
subplot(2,1,2);
plot(time, Freq2);
title('Area 2 Frequency Change vs Time')
xlabel('Time(s)')
ylabel('Frequency Change (pu)')
figure(5);
plot(time, TieLine);
title('Tie Line Power vs Time')
xlabel('Time(s)')
ylabel('Tie Line Power (pu)')

BIN
matlab_implementation/ace1.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 29 KiB

BIN
matlab_implementation/ace2.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 29 KiB

BIN
matlab_implementation/freq1.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 34 KiB

BIN
matlab_implementation/freq1_comparison.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 32 KiB

BIN
matlab_implementation/freq2.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 32 KiB

271
matlab_implementation/main.asv Normal file
View File

@@ -0,0 +1,271 @@
% Load disturbance for Area 1
PL1 = 0;
% Load disturbance for Area 2
PL2 = 0;
% Importing RBF Class
import RBF.*
% Defining an object (nn1) of the class RBF
nn1 = RBF( 0.13,0.21,0.25,0.9,0.98, 3);
% Defining an object (nn2) of the class RBF
nn2 = RBF( 0.13,0.21,0.25,0.9,0.98, 3);
% Defining and initializing variables
ut1 = []; % Variable for control signal for area 1
ut2 = []; % Variable for control signal for area 2
pl1 = []; % Variable for load disturbance for area 1
pl2 = []; % Variable for load disturbance for area 2
ACE1 = []; % Variable for area control error for area 1
ACE2 = []; % Variable for area control error for area 2
TieLine = []; % Variable for TieLine Power
Freq1 =[]; % Variable for Frequency change for area 2
Freq2 = [];% Variable for Frequency change for area 2
time = []; % Variable for time
K = [];
Tg = 0.08; % Initializing variable for governor time constant
Tp = 20; % Initializing variable for plant time constant
Tt = 0.3; % Initializing variable for Turbine time constant
Kp = 120; % Initializing variable for PID parameter Kp
T12 = 0.545/(2*pi); % Initializing variable for Tie line power constant
a12 = -1; % Initializing variable a12
R = 5; % Initializing variable for Regulation
T = 1/80; % Initializing Time step
dt = 1/80; % Initializing Time step
beta1 = 0.425; % Initializing variable for Frequency bias constant of area 1
beta2 = 0.425; % Initializing variable for Frequency bias constant of area 2
% Defining object System for the class TwoAreaPS
System = TwoAreaPS( Tg, Tp, Tt, Kp, T12, a12, R, T, beta1, beta2, [0,0,0],[0,0,0] );
% Initializing Ki
Ki = 0 ;
% Initializing Kd
Kd = 0 ;
% Initializing Kp
Kp = 0 ;
% Initializing weights
nn1.initialize_weights(1);
nn2.initialize_weights(2);
% Initializing control signals
ut_1 = 0;
ut_2 = 0;
% Duration of run
t = 200;
% Output variable
y = [];
for i = 1:t/dt
% Finding error signal e_t for area 1
e_t_a1 = 0 - System.yt_a1(1);
% Finding derivative of the output (delta y) for area 1
del_y_a1 = System.yt_a1(1) - System.yt_a1(2);
% Finding second derivate of the output (delta square y) for area 1
del2_y_a1 = System.yt_a1(1) - 2*System.yt_a1(2) + System.yt_a1(3);
% Defining the state X for area 1 controller
nn1.X = [ e_t_a1 ; -del_y_a1 ; -del2_y_a1];
% Calling HiddenLayer function using object nn1
nn1.HiddenLayer();
% Calling OutputLayer function using object nn1
nn1.OutputLayer();
% Finding error signal e_t for area 2
e_t_a2 = 0 - System.yt_a2(1);
% Finding derivative of the output (delta y) for area 2
del_y_a2 = System.yt_a2(1) - System.yt_a2(2);
% Finding second derivate of the output (delta square y) for area 2
del2_y_a2 = System.yt_a2(1) - 2*System.yt_a2(2) + System.yt_a2(3);
% Defining the state X for area 2 controller
nn2.X = [ e_t_a2 ; -del_y_a2 ; -del2_y_a2];
% Calling HiddenLayer function using object nn2
nn2.HiddenLayer();
% Calling OutputLayer function using object nn2
nn2.OutputLayer();
% Control signal (Equation 22 in the paper) for area 1
ut_1 = ut_1 + (nn1.K(2)*e_t_a1 - nn1.K(1)*del_y_a1 - nn1.K(3)*del2_y_a1);
% Control signal (Equation 22 in the paper) for area 2
ut_2 = ut_2 + (nn2.K(2)*e_t_a2 - nn2.K(1)*del_y_a2 - nn2.K(3)*del2_y_a2);
% Calculating delU1
delU1 = (nn1.K(2)*e_t_a1 - nn1.K(1)*del_y_a1 - nn1.K(3)*del2_y_a1);
K = [K nn1.K(3)];
% Calculating delU2
delU2 = (nn2.K(2)*e_t_a2 - nn2.K(1)*del_y_a2 - nn2.K(3)*del2_y_a2);
% Defining control signal ut1 for area 1
ut1 = [ut1 ut_1];
% Defining control signal ut2 for area 2
ut2 = [ut2 ut_2];
% Load disturbances
if (i*dt >=40 && i*dt <= 80)
PL1 = 0.003;
elseif( i*dt > 80 && i*dt <= 120)
PL1 = 0.006;
elseif( i*dt > 120 && i*dt <= 160)
PL1 = 0.009;
else
PL1 = 0;
end
PL2 = 0;
pl1 = [pl1 PL1];
pl2 = [pl2 PL2];
% Defining control signal Ut
Ut = [ ut_1 ; ut_2 ; PL1 ; PL2 ];
System.Output(Ut);
% Calling Update function using the object nn1 for area 1
nn1.Update(0 ,System.Y(1), System.yt_a1(1), System.yt_a1(2), System.yt_a1(3),delU1+1);
% Calling Update function using the object nn2 for area 2
nn2.Update(0 ,System.Y(2), System.yt_a2(1), System.yt_a2(2), System.yt_a2(3),delU2+1 );
% Defining ACE for area 1
ACE1 = [ACE1 System.Y(1)];
% Defining ACE for area 2
ACE2 = [ACE2 System.Y(2)];
% Defining Tie Line Power
TieLine = [ TieLine System.Y(3)];
% Defining Frequency Change for area 1
Freq1 = [ Freq1 System.Y(4)];
% Defining Frequency Change for area 2
Freq2 = [ Freq2 System.Y(5)];
% Defining variable time for our time period
time = [time i*dt];
% Metrics calculations
%Integral time absolute error(eq 38 in paper)
ITAE = [];
ITAE = [ITAE dt*(i*dt)*( abs(System.Y(1)) + abs(System.Y(2)))];
% Integral time squared error(eq 36 in paper)
ITSE = [];
ITSE = [ITSE dt*(i*dt)*( System.Y(1)^2 + System.Y(2)^2 )];
% Integral absolute error(eq 37 in paper)
IAE = [];
IAE = [IAE dt*( abs(System.Y(1)) + abs(System.Y(2)) ) ];
% Integral Square error(eq 39 in paper)
ISE = [];
ISE = [ISE dt*( System.Y(1)^2 + System.Y(2)^2 ) ];
end
% Displaying metrics values
fprintf('IAE is');
disp(sum(IAE));
fprintf('ISE is');
disp(sum(ISE));
fprintf('ITSE is');
disp(sum(ITSE));
fprintf('ITAE is');
disp(sum(ITAE));
% disp(K);
writematrix(ACE1,"ACE1.txt");
writematrix(ACE2,"ACE2.txt");
writematrix(time,"time.txt");
% plotting load graph
figure(1);
plot(time,pl1 ,time,pl2);
title(' Area 1 and Area 2 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
legend('Area 1', 'Area 2')
% ylabel('Load (pu)')
% subplot(2,1,2);
% plot(time, pl2);
% title('Area 2 Load vs Time')
% xlabel('Time(s)')
% plotting control signal graph
figure(2);
subplot(2,1,1);
plot(time,ut1 );
title(' Area 1 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
subplot(2,1,2);
plot(time, ut2);
title('Area 2 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
% plotting ACE graph
figure(3);
subplot(2,1,1);
plot(time, ACE1);
title(' Area 1 ACE vs Time for Adaptive PID')
xlabel('Time(s)')
ylabel('ACE (pu)')
subplot(2,1,2);
plot(time, ACE2);
title('Area 2 ACE vs Time Adaptive PID')
xlabel('Time(s)')
ylabel('ACE (pu)')
figure(4);
plot(time,K);
title(' Area 1 and Area 2 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
% legend('Area 1', 'Area 2')
figure(4);
subplot(2,1,1);
plot(time, Freq1);
title(' Area 1 Frequency Change vs Time')
xlabel('Time(s)')
ylabel('Frequency Change (pu)')
subplot(2,1,2);
plot(time, Freq2);
title('Area 2 Frequency Change vs Time')
xlabel('Time(s)')
ylabel('Frequency Change (pu)')
figure(5);
plot(time, TieLine);
title('Tie Line Power vs Time')
xlabel('Time(s)')
ylabel('Tie Line Power (pu)')

271
matlab_implementation/main.m Normal file
View File

@@ -0,0 +1,271 @@
% Load disturbance for Area 1
PL1 = 0;
% Load disturbance for Area 2
PL2 = 0;
% Importing RBF Class
import RBF.*
% Defining an object (nn1) of the class RBF
nn1 = RBF( 0.13,0.21,0.25,0.9,0.98, 3);
% Defining an object (nn2) of the class RBF
nn2 = RBF( 0.13,0.21,0.25,0.9,0.98, 3);
% Defining and initializing variables
ut1 = []; % Variable for control signal for area 1
ut2 = []; % Variable for control signal for area 2
pl1 = []; % Variable for load disturbance for area 1
pl2 = []; % Variable for load disturbance for area 2
ACE1 = []; % Variable for area control error for area 1
ACE2 = []; % Variable for area control error for area 2
TieLine = []; % Variable for TieLine Power
Freq1 =[]; % Variable for Frequency change for area 2
Freq2 = [];% Variable for Frequency change for area 2
time = []; % Variable for time
K = [];
Tg = 0.08; % Initializing variable for governor time constant
Tp = 20; % Initializing variable for plant time constant
Tt = 0.3; % Initializing variable for Turbine time constant
Kp = 120; % Initializing variable for PID parameter Kp
T12 = 0.545/(2*pi); % Initializing variable for Tie line power constant
a12 = -1; % Initializing variable a12
R = 5; % Initializing variable for Regulation
T = 1/80; % Initializing Time step
dt = 1/80; % Initializing Time step
beta1 = 0.425; % Initializing variable for Frequency bias constant of area 1
beta2 = 0.425; % Initializing variable for Frequency bias constant of area 2
% Defining object System for the class TwoAreaPS
System = TwoAreaPS( Tg, Tp, Tt, Kp, T12, a12, R, T, beta1, beta2, [0,0,0],[0,0,0] );
% Initializing Ki
Ki = 0 ;
% Initializing Kd
Kd = 0 ;
% Initializing Kp
Kp = 0 ;
% Initializing weights
nn1.initialize_weights(1);
nn2.initialize_weights(2);
% Initializing control signals
ut_1 = 0;
ut_2 = 0;
% Duration of run
t = 200;
% Output variable
y = [];
for i = 1:t/dt
% Finding error signal e_t for area 1
e_t_a1 = 0 - System.yt_a1(1);
% Finding derivative of the output (delta y) for area 1
del_y_a1 = System.yt_a1(1) - System.yt_a1(2);
% Finding second derivate of the output (delta square y) for area 1
del2_y_a1 = System.yt_a1(1) - 2*System.yt_a1(2) + System.yt_a1(3);
% Defining the state X for area 1 controller
nn1.X = [ e_t_a1 ; -del_y_a1 ; -del2_y_a1];
% Calling HiddenLayer function using object nn1
nn1.HiddenLayer();
% Calling OutputLayer function using object nn1
nn1.OutputLayer();
% Finding error signal e_t for area 2
e_t_a2 = 0 - System.yt_a2(1);
% Finding derivative of the output (delta y) for area 2
del_y_a2 = System.yt_a2(1) - System.yt_a2(2);
% Finding second derivate of the output (delta square y) for area 2
del2_y_a2 = System.yt_a2(1) - 2*System.yt_a2(2) + System.yt_a2(3);
% Defining the state X for area 2 controller
nn2.X = [ e_t_a2 ; -del_y_a2 ; -del2_y_a2];
% Calling HiddenLayer function using object nn2
nn2.HiddenLayer();
% Calling OutputLayer function using object nn2
nn2.OutputLayer();
% Control signal (Equation 22 in the paper) for area 1
ut_1 = ut_1 + (nn1.K(2)*e_t_a1 - nn1.K(1)*del_y_a1 - nn1.K(3)*del2_y_a1);
% Control signal (Equation 22 in the paper) for area 2
ut_2 = ut_2 + (nn2.K(2)*e_t_a2 - nn2.K(1)*del_y_a2 - nn2.K(3)*del2_y_a2);
% Calculating delU1
delU1 = (nn1.K(2)*e_t_a1 - nn1.K(1)*del_y_a1 - nn1.K(3)*del2_y_a1);
K = [K nn1.K(3)];
% Calculating delU2
delU2 = (nn2.K(2)*e_t_a2 - nn2.K(1)*del_y_a2 - nn2.K(3)*del2_y_a2);
% Defining control signal ut1 for area 1
ut1 = [ut1 ut_1];
% Defining control signal ut2 for area 2
ut2 = [ut2 ut_2];
% Load disturbances
if (i*dt >=40 && i*dt <= 80)
PL1 = 0.003;
elseif( i*dt > 80 && i*dt <= 120)
PL1 = 0.006;
elseif( i*dt > 120 && i*dt <= 160)
PL1 = 0.009;
else
PL1 = 0;
end
PL2 = 0;
pl1 = [pl1 PL1];
pl2 = [pl2 PL2];
% Defining control signal Ut
Ut = [ ut_1 ; ut_2 ; PL1 ; PL2 ];
System.Output(Ut);
% Calling Update function using the object nn1 for area 1
nn1.Update(0 ,System.Y(1), System.yt_a1(1), System.yt_a1(2), System.yt_a1(3),delU1+1);
% Calling Update function using the object nn2 for area 2
nn2.Update(0 ,System.Y(2), System.yt_a2(1), System.yt_a2(2), System.yt_a2(3),delU2+1 );
% Defining ACE for area 1
ACE1 = [ACE1 System.Y(1)];
% Defining ACE for area 2
ACE2 = [ACE2 System.Y(2)];
% Defining Tie Line Power
TieLine = [ TieLine System.Y(3)];
% Defining Frequency Change for area 1
Freq1 = [ Freq1 System.Y(4)];
% Defining Frequency Change for area 2
Freq2 = [ Freq2 System.Y(5)];
% Defining variable time for our time period
time = [time i*dt];
% Metrics calculations
%Integral time absolute error(eq 38 in paper)
ITAE = [];
ITAE = [ITAE dt*(i*dt)*( abs(System.Y(1)) + abs(System.Y(2)))];
% Integral time squared error(eq 36 in paper)
ITSE = [];
ITSE = [ITSE dt*(i*dt)*( System.Y(1)^2 + System.Y(2)^2 )];
% Integral absolute error(eq 37 in paper)
IAE = [];
IAE = [IAE dt*( abs(System.Y(1)) + abs(System.Y(2)) ) ];
% Integral Square error(eq 39 in paper)
ISE = [];
ISE = [ISE dt*( System.Y(1)^2 + System.Y(2)^2 ) ];
end
% Displaying metrics values
fprintf('IAE is');
disp(sum(IAE));
fprintf('ISE is');
disp(sum(ISE));
fprintf('ITSE is');
disp(sum(ITSE));
fprintf('ITAE is');
disp(sum(ITAE));
% disp(K);
writematrix(ACE1,"ACE1.txt");
writematrix(ACE2,"ACE2.txt");
writematrix(time,"time.txt");
% plotting load graph
figure(1);
plot(time,pl1 ,time,pl2);
title(' Area 1 and Area 2 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
legend('Area 1', 'Area 2')
% ylabel('Load (pu)')
% subplot(2,1,2);
% plot(time, pl2);
% title('Area 2 Load vs Time')
% xlabel('Time(s)')
% plotting control signal graph
figure(2);
subplot(2,1,1);
plot(time,ut1 );
title(' Area 1 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
subplot(2,1,2);
plot(time, ut2);
title('Area 2 Control Signal vs Time')
xlabel('Time(s)')
ylabel('Control Signal (pu)')
% plotting ACE graph
figure(3);
subplot(2,1,1);
plot(time, ACE1);
title(' Area 1 ACE vs Time for Adaptive PID')
xlabel('Time(s)')
ylabel('ACE (pu)')
subplot(2,1,2);
plot(time, ACE2);
title('Area 2 ACE vs Time Adaptive PID')
xlabel('Time(s)')
ylabel('ACE (pu)')
figure(4);
plot(time,K);
title(' Area 1 and Area 2 Load vs Time')
xlabel('Time(s)')
ylabel('Load (pu)')
% legend('Area 1', 'Area 2')
figure(4);
subplot(2,1,1);
plot(time, Freq1);
title(' Area 1 Frequency Change vs Time')
xlabel('Time(s)')
ylabel('Frequency Change (pu)')
subplot(2,1,2);
plot(time, Freq2);
title('Area 2 Frequency Change vs Time')
xlabel('Time(s)')
ylabel('Frequency Change (pu)')
figure(5);
plot(time, TieLine);
title('Tie Line Power vs Time')
xlabel('Time(s)')
ylabel('Tie Line Power (pu)')

3
matlab_implementation/mu_matrix1.txt Normal file
View File

@@ -0,0 +1,3 @@
1.845599726003620614e-02 1.516044794402831306e-03 8.995979242641155338e+00
5.780865431825088323e-03 4.713686580375981566e-03 3.674395894835833065e+00
9.482924956430698771e+00 -2.785259136223135812e+00 4.344232482720499000e+00

3
matlab_implementation/mu_matrix2.txt Normal file
View File

@@ -0,0 +1,3 @@
1.181781558224546773e-02 1.130746439805960235e-02 4.709530595006598652e+00
6.339533643002934378e-03 3.406172729208175026e-03 3.687143721330234047e+00
1.537852532877977518e+00 1.835690290426149396e+00 7.904398565729198189e+00

4
matlab_implementation/nonlinearmain.asv Normal file
View File

@@ -0,0 +1,4 @@
yd = [ 2.1*ones(1,100) 3.5*ones(1,100) 2*ones(1,100) 3*ones(1,100)];
[y,x,i] = nonlinearsystem(yd);
plot(x,y)

4
matlab_implementation/nonlinearmain.m Normal file
View File

@@ -0,0 +1,4 @@
yd = [ 2.1*ones(1,100) 3.5*ones(1,100) 2*ones(1,100) 3*ones(1,100)];
[y,x,i] = nonlinearsystem(yd);
plot(x,y)

50
matlab_implementation/nonlinearsystem.asv Normal file
View File

@@ -0,0 +1,50 @@
function [y, x, initial_states] = nonlinearsystem(yd)
import RBF.*
rbf = RBF(0.13,0.21,0.25,0.9,0.98,3);
t = 1;
yt_1 = 0.4;
yt_2 = 0.46;
yt_3 = 0.42;
initial_states = [yt_1, yt_2, yt_3];
dt = 1/400;
Ki = -0.07709546;
Kd = 0.58844546;
Kp = -0.01747239;
ut_1 = 0;
y = [];
x = [];
for i = 1:t/dt
e_t = yd(i) - yt_1;
del_y = yt_1 - yt_2;
del2_y = yt_1 - 2*yt_2 + yt_3;
rbf.X = [ e_t ; -del_y ; -del2_y];
rbf.HiddenLayer();
rbf.OutputLayer();
ut_1 = ut_1 + rbf.K(2)*e_t - rbf.K(1)*del_y - rbf.K(3)*del2_y;
disp(rbf.K);
% [y1,y2,y3] = deal(yt_1 , yt_2, yt_3);
[yt_1 , yt_2, yt_3] = deal(yt_1*yt_2*(yt_1 + 2.5) / ( 1 + yt_1^2 + yt_2^2 )+ ut_1 , yt_1 , yt_2) ;
y0 = yt_1;
rbf.Update(yd(i),y0);
y = [y, yt_1];
x = [x, i];
end
end

50
matlab_implementation/nonlinearsystem.m Normal file
View File

@@ -0,0 +1,50 @@
function [y, x, initial_states] = nonlinearsystem(yd)
import RBF.*
rbf = RBF(0.13,0.21,0.25,0.9,0.98,3);
t = 1;
yt_1 = 0.4;
yt_2 = 0.46;
yt_3 = 0.42;
initial_states = [yt_1, yt_2, yt_3];
dt = 1/400;
Ki = -0.07709546;
Kd = 0.58844546;
Kp = -0.01747239;
ut_1 = 0;
y = [];
x = [];
for i = 1:t/dt
e_t = yd(i) - yt_1;
del_y = yt_1 - yt_2;
del2_y = yt_1 - 2*yt_2 + yt_3;
rbf.X = [ e_t ; -del_y ; -del2_y];
rbf.HiddenLayer();
rbf.OutputLayer();
ut_1 = ut_1 + rbf.K(2)*e_t - rbf.K(1)*del_y - rbf.K(3)*del2_y;
disp(rbf.K);
% [y1,y2,y3] = deal(yt_1 , yt_2, yt_3);
[yt_1 , yt_2, yt_3] = deal(yt_1*yt_2*(yt_1 + 2.5) / ( 1 + yt_1^2 + yt_2^2 )+ ut_1 , yt_1 , yt_2) ;
y0 = yt_1;
rbf.Update(yd(i),y0);
y = [y, yt_1];
x = [x, i];
end
end

39
matlab_implementation/pso_with_two_area/Particle.asv Normal file
View File

@@ -0,0 +1,39 @@
classdef Particle
properties
fitness;
dim;
minx;
maxx;
position;
velocity;
best_part_pos;
best_part_fitnessVal;
rnd;
end
methods
function self = Particle(dim, minx, maxx)
rng('default')
self.rnd = rand();
self.position = zeros(1,dim);
self.velocity = zeros(1,dim);
self.best_part_pos = zeros(1,dim);
for i = 1:dim
self.position(i) = ((maxx - minx) * self.rnd + minx);
self.velocity(i) = ((maxx - minx) * self.rnd + minx);
self.fitness = system_with_pid(self.position) ;
self.best_part_pos = self.position;
self.best_part_fitnessVal = self.fitness ;
end
end
end
end

39
matlab_implementation/pso_with_two_area/Particle.m Normal file
View File

@@ -0,0 +1,39 @@
classdef Particle
properties
fitness;
dim;
minx;
maxx;
position;
velocity;
best_part_pos;
best_part_fitnessVal;
rnd;
end
methods
function self = Particle(dim, minx, maxx)
rng('default')
self.rnd = rand();
self.position = zeros(1,dim);
self.velocity = zeros(1,dim);
self.best_part_pos = zeros(1,dim);
for i = 1:dim
self.position(i) = ((maxx - minx) * self.rnd + minx);
self.velocity(i) = ((maxx - minx) * self.rnd + minx);
self.fitness = system_with_pid(self.position) ;
self.best_part_pos = self.position;
self.best_part_fitnessVal = self.fitness ;
end
end
end
end

65
matlab_implementation/pso_with_two_area/pso.asv Normal file
View File

@@ -0,0 +1,65 @@
function best_swarm_pos = pso(max_iter, n, dim, minx, maxx)
w = 0.729;
c1 = 1.49445;
c2 = 1.49445;
rng('default')
rnd = rand();
swarm = [];
for i = 1:n
swarm = [swarm, Particle(dim, minx, maxx)];
end
best_swarm_fitnessVal = Inf;
for i = 1:n
if swarm(i).fitness < best_swarm_fitnessVal
best_swarm_fitnessVal = swarm(i).fitness;
best_swarm_pos = swarm(i).position;
end
end
Iter = 0;
while Iter < max_iter
for i = 1:n
for k = 1:dim
r1 = rand();
r2 = rand();
swarm(i).velocity(k) = ((w * swarm(i).velocity(k)) + (c1 * r1 * (swarm(i).best_part_pos(k) - swarm(i).position(k))) + (c2 * r2 * (best_swarm_pos(k) -swarm(i).position(k)))) ;
if swarm(i).velocity(k) < minx
swarm(i).velocity(k) = minx;
elseif swarm(i).velocity(k) > maxx
swarm(i).velocity(k) = maxx;
end
for k = 1:dim
swarm(i).position(k) = swarm(i).position(k) + swarm(i).velocity(k);
end
swarm(i).fitness = system_with_pid(swarm(i).position);
if swarm(i).fitness < swarm(i).best_part_fitnessVal
swarm(i).best_part_fitnessVal = swarm(i).fitness;
swarm(i).best_part_pos = swarm(i).position;
end
if swarm(i).fitness < best_swarm_fitnessVal
best_swarm_fitnessVal = swarm(i).fitness;
best_swarm_pos = swarm(i).position;
end
end
end
Iter = Iter+1;
end
end

65
matlab_implementation/pso_with_two_area/pso.m Normal file
View File

@@ -0,0 +1,65 @@
function best_swarm_pos = pso(max_iter, n, dim, minx, maxx)
w = 0.729;
c1 = 1.49445;
c2 = 1.49445;
rng('default')
rnd = rand();
swarm = [];
for i = 1:n
swarm = [swarm, Particle(dim, minx, maxx)];
end
best_swarm_fitnessVal = Inf;
for i = 1:n
if swarm(i).fitness < best_swarm_fitnessVal
best_swarm_fitnessVal = swarm(i).fitness;
best_swarm_pos = swarm(i).position;
end
end
Iter = 0;
while Iter < max_iter
for i = 1:n
for k = 1:dim
r1 = rand();
r2 = rand();
swarm(i).velocity(k) = ((w * swarm(i).velocity(k)) + (c1 * r1 * (swarm(i).best_part_pos(k) - swarm(i).position(k))) + (c2 * r2 * (best_swarm_pos(k) -swarm(i).position(k)))) ;
if swarm(i).velocity(k) < minx
swarm(i).velocity(k) = minx;
elseif swarm(i).velocity(k) > maxx
swarm(i).velocity(k) = maxx;
end
for k = 1:dim
swarm(i).position(k) = swarm(i).position(k) + swarm(i).velocity(k);
end
swarm(i).fitness = system_with_pid(swarm(i).position);
if swarm(i).fitness < swarm(i).best_part_fitnessVal
swarm(i).best_part_fitnessVal = swarm(i).fitness;
swarm(i).best_part_pos = swarm(i).position;
end
if swarm(i).fitness < best_swarm_fitnessVal
best_swarm_fitnessVal = swarm(i).fitness;
best_swarm_pos = swarm(i).position;
end
end
end
Iter = Iter+1;
end
end

63
matlab_implementation/pso_with_two_area/pso_example.asv Normal file
View File

@@ -0,0 +1,63 @@
% Initialization
% Parameters
clear
clc
iterations = 50;
W = 0.9;
C1 = 2.0;
C2 = 2.0;
n = 49;
% ---- initial swarm position -----
index = 1;
for i = 1 : 7
for j = 1 : 7
particle(index, 1, 1) = i;
particle(index, 1, 2) = j;
index = index + 1;
end
end
particle(:, 4, 1) = 1000; % best value so far
particle(:, 2, :) = 0; % initial velocity
%% Iterations
for iter = 1 : iterations
%-- evaluating position & quality ---
for i = 1 : n
particle(i, 1, 1) = particle(i, 1, 1) + particle(i, 2, 1)/1.3; %update x position
particle(i, 1, 2) = particle(i, 1, 2) + particle(i, 2, 2)/1.3; %update y position
Kp = particle(i, 1, 1);
Ki = particle(i, 1, 2);
Gc=tf([Kp Ki],[1 0]);
G=tf(1,[1 2 3]);
Q=feedback(Gc*G,1);
Z=stepinfo(Q);
a=Z.RiseTime;
b=Z.SettlingTime;
% Z.Overshoot=350;
c=Z.Overshoot;
d=Z.Undershoot;
% e=Z.SettlingMax;
% alpha=0.5;
F = min((a)+min(c)+min(b)+min(d)); % fitness evaluation
if F < particle(i, 4, 1) % if new position is better
particle(i, 3, 1) = particle(i, 1, 1); % update best x,
particle(i, 3, 2) = particle(i, 1, 2); % best y postions
particle(i, 4, 1) = F; % and best value
end
end
[temp, gbest] = min(particle(:, 4, 1)); % global best position
%--- updating velocity vectors
for i = 1 : n
particle(i, 2, 1) = rand*W*particle(i, 2, 1) + C1*rand*(particle(i, 3, 1) - particle(i, 1, 1)) + C2*rand*(particle(gbest, 3, 1) - particle(i, 1, 1)); %x velocity component
particle(i, 2, 2) = rand*W*particle(i, 2, 2) + C1*rand*(particle(i, 3, 2) - particle(i, 1, 2)) + C2*rand*(particle(gbest, 3, 2) - particle(i, 1, 2)); %y velocity component
end
%% Plotting the swarm
disp(gbest)
clf
plot(particle(:, 1, 1), particle(:, 1, 2), 'x')
axis([-2 10 -2 10]);
pause(0)
end
disp('best value is')
disp(gbest)

63
matlab_implementation/pso_with_two_area/pso_example.m Normal file
View File

@@ -0,0 +1,63 @@
% Initialization
% Parameters
clear
clc
iterations = 50;
W = 0.9;
C1 = 2.0;
C2 = 2.0;
n = 49;
% ---- initial swarm position -----
index = 1;
for i = 1 : 7
for j = 1 : 7
particle(index, 1, 1) = i;
particle(index, 1, 2) = j;
index = index + 1;
end
end
particle(:, 4, 1) = 1000; % best value so far
particle(:, 2, :) = 0; % initial velocity
%% Iterations
for iter = 1 : iterations
%-- evaluating position & quality ---
for i = 1 : n
particle(i, 1, 1) = particle(i, 1, 1) + particle(i, 2, 1)/1.3; %update x position
particle(i, 1, 2) = particle(i, 1, 2) + particle(i, 2, 2)/1.3; %update y position
Kp = particle(i, 1, 1);
Ki = particle(i, 1, 2);
Gc=tf([Kp Ki],[1 0]);
G=tf(1,[1 2 3]);
Q=feedback(Gc*G,1);
Z=stepinfo(Q);
a=Z.RiseTime;
b=Z.SettlingTime;
% Z.Overshoot=350;
c=Z.Overshoot;
d=Z.Undershoot;
% e=Z.SettlingMax;
% alpha=0.5;
F = min((a)+min(c)+min(b)+min(d)); % fitness evaluation
if F < particle(i, 4, 1) % if new position is better
particle(i, 3, 1) = particle(i, 1, 1); % update best x,
particle(i, 3, 2) = particle(i, 1, 2); % best y postions
particle(i, 4, 1) = F; % and best value
end
end
[temp, gbest] = min(particle(:, 4, 1)); % global best position
%--- updating velocity vectors
for i = 1 : n
particle(i, 2, 1) = rand*W*particle(i, 2, 1) + C1*rand*(particle(i, 3, 1) - particle(i, 1, 1)) + C2*rand*(particle(gbest, 3, 1) - particle(i, 1, 1)); %x velocity component
particle(i, 2, 2) = rand*W*particle(i, 2, 2) + C1*rand*(particle(i, 3, 2) - particle(i, 1, 2)) + C2*rand*(particle(gbest, 3, 2) - particle(i, 1, 2)); %y velocity component
end
%% Plotting the swarm
disp(gbest)
clf
plot(particle(:, 1, 1), particle(:, 1, 2), 'x')
axis([-2 10 -2 10]);
pause(0)
end
disp('best value is')
disp(gbest)

114
matlab_implementation/pso_with_two_area/pso_example_2.asv Normal file
View File

@@ -0,0 +1,114 @@
clear
clc
num_particles = 50; % Size of the swarm " no of birds "
max_iter = 50; % Maximum number of "birds steps"
dim = 6; % Dimension of the problem
c2 = 1.49445; % PSO parameter C1
c1 = 1.49445; % PSO parameter C2
w = 0.729; % pso momentum or inertia
fitness = zeros(num_particles,max_iter);
%-----------------------------%
% initialize the parameter %
%-----------------------------%
R1 = rand(dim, num_particles);
R2 = rand(dim, num_particles);
current_fitness = zeros(num_particles,1);
%------------------------------------------------%
% Initializing swarm and velocities and position %
%------------------------------------------------%
current_position = 10 * (rand(dim, num_particles)-0.5);
velocity = 0.3 * randn(dim, num_particles) ;
local_best_position = current_position ;
%-------------------------------------------%
% Evaluate initial population %
%-------------------------------------------%
for i = 1 : num_particles
current_fitness(i) = system_with_pid(current_position(:,i));
end
local_best_fitness = current_fitness ;
[global_best_fitness,g] = min(local_best_fitness) ;
for i = 1 : num_particles
global_best_position(:,i) = local_best_position(:,g) ;
end
% %-------------------%
% % VELOCITY UPDATE %
% %-------------------%
%
% velocity = w *velocity + c1*(R1.*(local_best_position-current_position)) + c2*(R2.*(global_best_position-current_position));
%
% %------------------%
% % SWARMUPDATE %
% %------------------%
%
%
% current_position = current_position + velocity ;
%------------------------%
% evaluate a new swarm %
%------------------------%
iter = 0 ; % Iterations’counter
while ( iter < max_iter )
if mod(iter,10)==0 && iter>1
% fprintf('Iter = %d', iter);
% fprintf('Best fitness = %.7f', current_global_best_fitness);
disp(iter);
disp(current_global_best_fitness);
end
iter = iter + 1;
for i = 1 : num_particles
current_fitness(i) = system_with_pid(current_position(:,i)) ;
end
for i = 1 : num_particles
if current_fitness(i) < local_best_fitness(i)
local_best_fitness(i) = current_fitness(i);
local_best_position(:,i) = current_position(:,i) ;
end
fitness(i,iter) = local_best_fitness(i);
end
[current_global_best_fitness,g] = min(local_best_fitness);
if current_global_best_fitness < global_best_fitness
global_best_fitness = current_global_best_fitness;
for i = 1 : num_particles
global_best_position(:,i) = local_best_position(:,g);
end
end
velocity = w *velocity + c1*(R1.*(local_best_position-current_position)) + c2*(R2.*(global_best_position-current_position));
current_position = current_position + velocity;
end % end of while loop its mean the end of all step that the birds move it
xx = fitness(:,max_iter);
% disp(current_global_best_fitness);
disp(xx);
[Y,I] = min(xx);
% disp(Y);
disp(current_position(:,I));
%

114
matlab_implementation/pso_with_two_area/pso_example_2.m Normal file
View File

@@ -0,0 +1,114 @@
clear
clc
num_particles = 50; % Size of the swarm " no of birds "
max_iter = 50; % Maximum number of "birds steps"
dim = 6; % Dimension of the problem
c2 = 1.49445; % PSO parameter C1
c1 = 1.49445; % PSO parameter C2
w = 0.729; % pso momentum or inertia
fitness = zeros(num_particles,max_iter);
%-----------------------------%
% initialize the parameter %
%-----------------------------%
R1 = rand(dim, num_particles);
R2 = rand(dim, num_particles);
current_fitness = zeros(num_particles,1);
%------------------------------------------------%
% Initializing swarm and velocities and position %
%------------------------------------------------%
current_position = 10 * (rand(dim, num_particles)-0.5);
velocity = 0.3 * randn(dim, num_particles) ;
local_best_position = current_position ;
%-------------------------------------------%
% Evaluate initial population %
%-------------------------------------------%
for i = 1 : num_particles
current_fitness(i) = system_with_pid(current_position(:,i));
end
local_best_fitness = current_fitness ;
[global_best_fitness,g] = min(local_best_fitness) ;
for i = 1 : num_particles
global_best_position(:,i) = local_best_position(:,g) ;
end
% %-------------------%
% % VELOCITY UPDATE %
% %-------------------%
%
% velocity = w *velocity + c1*(R1.*(local_best_position-current_position)) + c2*(R2.*(global_best_position-current_position));
%
% %------------------%
% % SWARMUPDATE %
% %------------------%
%
%
% current_position = current_position + velocity ;
%------------------------%
% evaluate a new swarm %
%------------------------%
iter = 0 ; % Iterations’counter
while ( iter < max_iter )
if mod(iter,10)==0 && iter>1
% fprintf('Iter = %d', iter);
% fprintf('Best fitness = %.7f', current_global_best_fitness);
disp(iter);
disp(current_global_best_fitness);
end
iter = iter + 1;
for i = 1 : num_particles
current_fitness(i) = system_with_pid(current_position(:,i)) ;
end
for i = 1 : num_particles
if current_fitness(i) < local_best_fitness(i)
local_best_fitness(i) = current_fitness(i);
local_best_position(:,i) = current_position(:,i) ;
end
fitness(i,iter) = local_best_fitness(i);
end
[current_global_best_fitness,g] = min(local_best_fitness);
if current_global_best_fitness < global_best_fitness
global_best_fitness = current_global_best_fitness;
for i = 1 : num_particles
global_best_position(:,i) = local_best_position(:,g);
end
end
velocity = w *velocity + c1*(R1.*(local_best_position-current_position)) + c2*(R2.*(global_best_position-current_position));
current_position = current_position + velocity;
end % end of while loop its mean the end of all step that the birds move it
xx = fitness(:,max_iter);
% disp(current_global_best_fitness);
disp(xx);
[Y,I] = min(xx);
% disp(Y);
disp(current_position(:,I));
%

12
matlab_implementation/pso_with_two_area/pso_main.asv Normal file
View File

@@ -0,0 +1,12 @@
dim = 6;
% fitness = @system_with_pid;
num_particles = 40;
max_iter = 50;
best_position = pso(max_iter, num_particles, dim, -1, 1);
disp(best_position);
fitnessVal = system_with_pid(best_position);
disp(fitnessVal);

11
matlab_implementation/pso_with_two_area/pso_main.m Normal file
View File

@@ -0,0 +1,11 @@
dim = 6;
num_particles = 40;
max_iter = 50;
best_position = pso(max_iter, num_particles, dim, -1, 1);
disp(best_position);
fitnessVal = system_with_pid(best_position);
disp(fitnessVal);

87
matlab_implementation/pso_with_two_area/system_with_pid.asv Normal file
View File

@@ -0,0 +1,87 @@
function fitness = system_with_pid(k_values)
PL1 = 0;
PL2 = 0;
import RBF.*
ut1 = [];
ut2 = [];
pl1 = [];
pl2 = [];
ACE1 = [];
ACE2 = [];
time = [];
Tg = 0.08;
Tp = 20;
Tt = 0.3;
Kp = 120;
T12 = 0.545/(2*pi);
a12 = -1;
R = 5;
T = 1/80;
dt = 1/80;
beta1 = 0.425;
beta2 = 0.425;
System = TwoAreaPS( Tg, Tp, Tt, Kp, T12, a12, R, T, beta1, beta2, [0,0,0],[0,0,0] );
KP1 = k_values(1);
KI1 = k_values(2);
KD1 = k_values(3);
KP2 = k_values(4);
KI2 = k_values(5);
KD2 = k_values(6);
ut_1 = 0;
ut_2 = 0;
t = 200;
y = [];
x = [];
for i = 1:t/dt
e_t_a1 = 0 - System.yt_a1(1);
del_y_a1 = System.yt_a1(1) - System.yt_a1(2);
del2_y_a1 = System.yt_a1(1) - 2*System.yt_a1(2) + System.yt_a1(3);
e_t_a2 = 0 - System.yt_a2(1);
del_y_a2 = System.yt_a2(1) - System.yt_a2(2);
del2_y_a2 = System.yt_a2(1) - 2*System.yt_a2(2) + System.yt_a2(3);
ut_1 = ut_1 + KI1*e_t_a1 - KP1*del_y_a1 - KD1*del2_y_a1;
ut_2 = ut_2 + KI2*e_t_a2 - KP2*del_y_a2 - KD2*del2_y_a2;
ut1 = [ut1 ut_1];
ut2 = [ut2 ut_2];
if (i*dt >=40 && i*dt <= 80)
PL1 = 0.003;
elseif( i*dt > 80 && i*dt <= 120)
PL1 = 0.006;
elseif( i*dt > 120 && i*dt <= 160)
PL1 = 0.009;
else
PL1 = 0;
end
PL2 = 0;
pl1 = [pl1 PL1];
pl2 = [pl2 PL2];
Ut = [ ut_1 ; ut_2 ; PL1 ; PL2 ];
System.Output(Ut);
performance_index = [];
performance_index = [performance_index, dt*(i*dt)*( abs(System.Y(1)) + abs(System.Y(2)))];
end
fitness = sum(performance_index);
end

87
matlab_implementation/pso_with_two_area/system_with_pid.m Normal file
View File

@@ -0,0 +1,87 @@
function fitness = system_with_pid(k_values)
PL1 = 0;
PL2 = 0;
import RBF.*
ut1 = [];
ut2 = [];
pl1 = [];
pl2 = [];
ACE1 = [];
ACE2 = [];
time = [];
Tg = 0.08;
Tp = 20;
Tt = 0.3;
Kp = 120;
T12 = 0.545/(2*pi);
a12 = -1;
R = 5;
T = 1/80;
dt = 1/80;
beta1 = 0.425;
beta2 = 0.425;
System = TwoAreaPS( Tg, Tp, Tt, Kp, T12, a12, R, T, beta1, beta2, [0,0,0],[0,0,0] );
KP1 = k_values(1);
KI1 = k_values(2);
KD1 = k_values(3);
KP2 = k_values(4);
KI2 = k_values(5);
KD2 = k_values(6);
ut_1 = 0;
ut_2 = 0;
t = 200;
y = [];
x = [];
for i = 1:t/dt
e_t_a1 = 0 - System.yt_a1(1);
del_y_a1 = System.yt_a1(1) - System.yt_a1(2);
del2_y_a1 = System.yt_a1(1) - 2*System.yt_a1(2) + System.yt_a1(3);
e_t_a2 = 0 - System.yt_a2(1);
del_y_a2 = System.yt_a2(1) - System.yt_a2(2);
del2_y_a2 = System.yt_a2(1) - 2*System.yt_a2(2) + System.yt_a2(3);
ut_1 = ut_1 + KI1*e_t_a1 - KP1*del_y_a1 - KD1*del2_y_a1;
ut_2 = ut_2 + KI2*e_t_a2 - KP2*del_y_a2 - KD2*del2_y_a2;
ut1 = [ut1 ut_1];
ut2 = [ut2 ut_2];
if (i*dt >=40 && i*dt <= 80)
PL1 = 0.003;
elseif( i*dt > 80 && i*dt <= 120)
PL1 = 0.006;
elseif( i*dt > 120 && i*dt <= 160)
PL1 = 0.009;
else
PL1 = 0;
end
PL2 = 0;
pl1 = [pl1 PL1];
pl2 = [pl2 PL2];
Ut = [ ut_1 ; ut_2 ; PL1 ; PL2 ];
System.Output(Ut);
performance_index = [];
performance_index = [performance_index, dt*(i*dt)*( abs(System.Y(1)) + abs(System.Y(2)))];
end
fitness = sum(performance_index);
end

1
matlab_implementation/sigma_matrix1.txt Normal file
View File

@@ -0,0 +1 @@
1.958288059466925191e+00 2.405280511179789561e-03 4.456775142277715318e+00

1
matlab_implementation/sigma_matrix2.txt Normal file
View File

@@ -0,0 +1 @@
4.844955477517413733e+00 1.878559667561860547e+00 4.493629460068961912e+00

BIN
matlab_implementation/tieline.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 33 KiB

1
matlab_implementation/time.txt Normal file
View File

File diff suppressed because one or more lines are too long

61
matlab_implementation/twoarea_pso.asv Normal file
View File

@@ -0,0 +1,61 @@
% Initialization
% Parameters
clear
clc
iterations = 50;
W = 0.9;
C1 = 2.0;
C2 = 2.0;
n = 49;
% ---- initial swarm position -----
index = 1;
for i = 1 : 7
for j = 1 : 7
particle(index, 1, 1) = i;
particle(index, 1, 2) = j;
index = index + 1;
end
end
particle(:, 4, 1) = 1000; % best value so far
particle(:, 2, :) = 0; % initial velocity
%% Iterations
for iter = 1 : iterations
%-- evaluating position & quality ---
for i = 1 : n
particle(i, 1, 1) = particle(i, 1, 1) + particle(i, 2, 1)/1.3; %update x position
particle(i, 1, 2) = particle(i, 1, 2) + particle(i, 2, 2)/1.3; %update y position
Kp = particle(i, 1, 1);
Ki = particle(i, 1, 2);
Gc=tf([Kp Ki],[1 0]);
G=tf(1,[1 2 3]);
Q=feedback(Gc*G,1);
Z=stepinfo(Q);
a=Z.RiseTime;
b=Z.SettlingTime;
% Z.Overshoot=350;
c=Z.Overshoot;
d=Z.Undershoot;
% e=Z.SettlingMax;
% alpha=0.5;
F = min((a)+min(c)+min(b)+min(d)); % fitness evaluation
if F < particle(i, 4, 1) % if new position is better
particle(i, 3, 1) = particle(i, 1, 1); % update best x,
particle(i, 3, 2) = particle(i, 1, 2); % best y postions
particle(i, 4, 1) = F; % and best value
end
end
[temp, gbest] = min(particle(:, 4, 1)); % global best position
%--- updating velocity vectors
for i = 1 : n
particle(i, 2, 1) = rand*W*particle(i, 2, 1) + C1*rand*(particle(i, 3, 1) - particle(i, 1, 1)) + C2*rand*(particle(gbest, 3, 1) - particle(i, 1, 1)); %x velocity component
particle(i, 2, 2) = rand*W*particle(i, 2, 2) + C1*rand*(particle(i, 3, 2) - particle(i, 1, 2)) + C2*rand*(particle(gbest, 3, 2) - particle(i, 1, 2)); %y velocity component
end
%% Plotting the swarm
disp(gbest)
clf
plot(particle(:, 1, 1), particle(:, 1, 2), 'x')
axis([-2 150 -2 150]);
pause(0)
end

BIN
matlab_implementation/untitled8.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 17 KiB

BIN
matlab_implementation/untitled9.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 17 KiB

Some files were not shown because too many files have changed in this diff Show More