Initial commit

LICENSE.txt

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2021 Georgios S. Misyris, Jochen Stiasny, Spyros Chatzivasileiadis
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

@@ -0,0 +1,32 @@
+# Capturing Power System Dynamics by Physics-Informed Neural Networks and Optimization
+
+This repository is the official implementation of [Capturing Power System Dynamics by Physics-Informed Neural Networks and Optimization](https://arxiv.org/abs/2103.17004). 
+
+## Environment
+
+To install and activate the environment using conda run:
+
+```setup
+conda env create -f environment.yml
+conda activate cdc_2021_converter_PINN
+```
+
+## Code structure
+
+The repository is devided into three main sections:
+- data_simulation
+- neural_network_training
+- neural_network_verification
+
+The paper is focused on the training and verification, the data simulation is added for completeness, the files data_simulation/training_data.mat and data_simulation/validation_data.mat stem from the simulation.
+
+## Neural network training
+The file training_process.py acts as the central file that combines the data preparation, the model setup, model training, and model saving, as well as a basic visualisation. The file model_weihgts.h5 contains a set of trained weights that can be in the model and the folder 'logs' shows the corresponding training process. 
+
+## Neural network verification
+The csv files contain the trained weights and biases that are subsequently used in the verification. In a first step 'Tighten_ReLU_Bounds_MILP.m' computes tighter bounds of the ReLU units which are then stored in 'zk_hat_max.mat' and 'zk_hat_min.mat'. The files starting with 'Compute_Maximum_' contain the optimization problems that are solved in order to analyse the system. The files 'system_analysis_power_loss_based.m' and 'system_analysis_voltage_deviation_based.m' utilize these to produce Fig. 4. and Fig. 5 in the paper.  
+
+## References
+The implementation of the Physics-Informed Neural Networks is done in [TensorFlow](https://www.tensorflow.org) (Martín Abadi et al., TensorFlow: Large-scale machine learning on heterogeneous systems, 2015. Software available from tensorflow.org.). 
+The neural network verification is an adaption from [Verification of neural networkbehaviour: Formal guarantees for power system applications](https://ieeexplore.ieee.org/abstract/document/9141308) (A. Venzke and S. Chatzivasileiadis, “Verification of neural networkbehaviour: Formal guarantees for power system applications,”IEEETransactions on Smart Grid, vol. 12, no. 1, pp. 383–397, 2021.)
+We furthermore use [Gurobi](https://www.gurobi.com) and [YALMIP](https://yalmip.github.io/) (J. Löfberg, “Yalmip : A toolbox for modeling and optimization in matlab,” in In Proceedings of the CACSD Conference, Taipei, Taiwan, 2004.) for solving the resulting optimization problems.

data_simulation/Power Flow Initialization/CheckTolerance.m

@@ -0,0 +1,14 @@
+function success = CheckTolerance(F,n)
+
+normF = norm(F, inf);
+
+global max_iter 
+global errortol
+
+if normF<errortol && n<max_iter
+   success=1;
+else
+   success=0;
+end
+
+end

data_simulation/Power Flow Initialization/DisplayResults.m

@@ -0,0 +1,39 @@
+function DisplayResults(V,Ybus,Y_from,Y_to,br_f,br_t,buscode)
+
+ n=length(buscode);
+
+    S_to   = V(br_t).*conj(Y_to*V);
+    S_from = V(br_f).*conj(Y_from*V); 
+
+
+S_inj = V.*conj(Ybus*V);
+
+
+ fprintf(' \n ============================================================================= \n'); 
+ fprintf('                            BUS RESULTS'); 
+ fprintf(' \n ============================================================================= \n');
+
+ fprintf(' \n Bus              Voltage                Generation               Load         \n');
+ fprintf(' \n  #          Mag(pu)  Ang(deg)         P(pu)    Q(pu)        P(pu)       Q(pu) \n')
+ fprintf(' \n ----        -----------------         ---------------       ----------------- \n');
+ for i=1:n
+ if real(S_inj(i))>0  
+ fprintf('\n  %d           %.3f     %.2f           %.2f       %.2f        -             -   \n', i , abs(V(i)), angle(V(i))*180/pi,  real(S_inj(i)) ,imag(S_inj(i)));
+ end
+ if real(S_inj(i))<=0
+ fprintf('\n  %d           %.3f     %.2f            -          -         %.2f        %.2f   \n', i , abs(V(i)), angle(V(i))*180/pi,  abs(real(S_inj(i))) ,(imag(S_inj(i))));
+ end
+ end
+
+ fprintf(' \n ============================================================================= \n'); 
+ fprintf('                                    BRANCH FLOW'); 
+ fprintf(' \n ============================================================================= \n');
+
+ fprintf(' \n Branch   From    To       From Bus   Injection        To Bus        Injection \n');
+ fprintf(' \n   #      Bus     Bus      P(pu)      Q(pu)            P(pu)         Q(pu)     \n')
+ fprintf(' \n ------  ------  -----     --------  ----------      --------       ---------- \n');
+ for i=1:n-1
+ fprintf(' \n   %d      %d     %.d           %.2f          %.2f        %.2f           %.2f\n',i,br_f(i),br_t(i),real(S_from(i)),imag(S_from(i)),real(S_to(i)),imag(S_to(i)));
+ end
+
+end

data_simulation/Power Flow Initialization/DisplayResults_2.m

@@ -0,0 +1,39 @@
+function DisplayResults(V,Ybus,Y_from,Y_to,br_f,br_t,buscode)
+
+ n=length(buscode);
+
+    S_to   = V(br_t).*conj(Y_to*V)
+    S_from = V(br_f).*conj(Y_from*V) 
+
+
+S_inj = V.*conj(Ybus*V)
+
+
+ fprintf(' \n ======================================================='); 
+ fprintf(' \n                    BUS RESULTS'); 
+ fprintf(' \n ======================================================= \n');
+
+ fprintf(' \n Bus       Voltage       Generation         Load\n');
+ fprintf(' \n  #   Mag(pu)  Ang(deg)  P(pu)  Q(pu)    P(pu)  Q(pu)\n')
+ fprintf(' \n ---- -----------------  ------------    ------------ \n');
+ for i=1:n
+ if real(S_inj(i))>0  
+ fprintf('   %d    %.3f     %.2f     %.2f   %.2f      -     -\n', i , abs(V(i)), angle(V(i))*180/pi,  real(S_inj(i)) ,imag(S_inj(i)));
+ end
+ if real(S_inj(i))<=0
+ fprintf('   %d    %.3f     %.2f      -      -      %.2f   %.2f\n', i , abs(V(i)), angle(V(i))*180/pi,  abs(real(S_inj(i))) ,abs(imag(S_inj(i))));
+ end
+ end
+
+ fprintf(' \n ====================================================================='); 
+ fprintf(' \n                            BRANCH FLOW'); 
+ fprintf(' \n ===================================================================== \n');
+
+ fprintf(' \n Branch   From    To      From Bus  Injection    To Bus    Injection\n');
+ fprintf(' \n   #      Bus     Bus      P(pu)      Q(pu)       P(pu)      Q(pu)\n')
+ fprintf(' \n ------  ------  -----    --------  ----------  --------   ----------\n');
+
+ fprintf(' \n  %d         %d       %.d      %.2f       %.2f        %.2f      %.2f\n',i,br_f,br_t,real(S_from),imag(S_from),real(S_to),imag(S_to));
+
+
+end

data_simulation/Power Flow Initialization/Plot_losses.m

@@ -0,0 +1,6 @@
+
+c = categorical({'0.2 f_{n}','0.5 f_{n}','0.8 f_{n}', 'f_{n}' ,' f_{n}'});
+PQlos = [Plosses(1) Qlosses(1);Plosses(2) Qlosses(2);Plosses(3) Qlosses(3);Plosses(4) Qlosses(4);Plosses(5) Qlosses(5)];
+
+figure()
+bar(n,PQlos)

data_simulation/Power Flow Initialization/PowerFlowNewton.m

@@ -0,0 +1,45 @@
+function [V,success,n]=PowerFlowNewton(Ybus,Sbus,V0,ref,pv_index,pq_index)
+
+global max_iter
+success=0;
+n=0;   %initialize a flag and a counter
+V=real(V0)+imag(V0)*1i; % Complex bus voltage vector
+% evaluate F(x0)
+
+F = calculate_F(Ybus,Sbus,V,pv_index,pq_index);
+
+fprintf('\n iteration   max P and Q mismatch (p.u.)');
+fprintf(' \n --------  ---------------------------');
+
+F
+%Check wether tolerance is reached
+success = CheckTolerance(F,n);
+
+
+% Start the Newton Iteration
+
+while (~success) && (n < max_iter)
+   n=n+1; % update iteration counter
+   % Generate the Jacobian
+   [J_dS_dVm, J_dS_dTheta] = generate_Derivatives(Ybus, V);
+   J = generate_Jacobian(J_dS_dVm, J_dS_dTheta, pv_index,pq_index);
+   % Compute the update step
+   dx = (J \ F);
+   %Update the voltages, calculate F and check whether tolerence is reached
+   [V,Vm,Theta] = Update_Voltages(dx,V,pv_index,pq_index);
+   F = calculate_F(Ybus,Sbus,V,pv_index,pq_index);
+   success = CheckTolerance(F,n); 
+
+end
+
+if ~success
+fprintf('\n Newton''s method did not converge in %d iterations.\n', n);
+end
+
+
+if success
+fprintf('\n OK.\n', n);
+end
+
+
+end

data_simulation/Power Flow Initialization/Power_flow_1VSC.m

@@ -0,0 +1,19 @@
+
+% clear all;
+% close all;
+% clc;
+
+global max_iter 
+global errortol
+max_iter=100000; 
+errortol=1e-3;
+
+% Load the Network Data, carry out the power flow and display results
+System_matrix_1;
+[V, success, n] = PowerFlowNewton(Ybus,Sbus,V0,ref,pv_index,pq_index);
+if (success)   
+    DisplayResults(V,Ybus,Y_from,Y_to,br_f,br_t,buscode);
+end
+S_inj = V.*conj(Ybus*V);
+    S_to   = V(br_t).*conj(Y_to*V);
+    S_from = V(br_f).*conj(Y_from*V); 

data_simulation/Power Flow Initialization/System_matrix_1.m

@@ -0,0 +1,37 @@
+
+Zfvsc1 = R+1i*L; Yfvsc1=1/Zfvsc1;
+
+
+
+
+
+Ybus = [Yfvsc1       -Yfvsc1 ;
+    -Yfvsc1          Yfvsc1];
+
+
+Sbus = [Pref+1i*Qref; 0+1i*0];
+V0 = 1*ones(length(Sbus),1);
+V0(1) = 1;
+V0(2) = 1; 
+buscode = [1; 3];
+
+pq_index = find(buscode==1);% Find indices for all PQ-busses
+pv_index = find(buscode==2); % Find indices for all PV-busses
+ref = find(buscode==3); % Find index for ref bus
+
+% Create branch matrices
+n_br = 1; n_bus = 2; % number of branches , number of busses
+Y_from = zeros(n_br,n_bus); % Create the two branch admittance matrices
+Y_to = zeros(n_br,n_bus);
+
+br_f = [1]; % The from busses
+br_t = [2]; % The to busses
+br_Y = [Yfvsc1]; % The series admittance of each branch
+
+for k = 1:length(br_f) % Fill in the matrices
+Y_from(k,br_f(k)) = br_Y(k);
+Y_from(k,br_t(k)) = -br_Y(k);
+Y_to(k,br_f(k)) = -br_Y(k);
+Y_to(k,br_t(k)) = br_Y(k);
+end
+

data_simulation/Power Flow Initialization/Update_Voltages.m

@@ -0,0 +1,24 @@
+function [V,Vm,Theta] = Update_Voltages(dx,V,pv_index,pq_index)
+
+N1 = 1; N2 = length(pv_index); %% dx(N1:N2)-ang. on the pv buses
+N3 = N2 + 1; N4 = N2+length(pq_index); %% dx(N3:N4)-ang. on the pq buses
+N5 = N4 + 1; N6 = N4+length(pq_index); %% dx(N5:N6)-mag. on the pq buses
+
+
+
+Theta = angle(V); 
+Vm = abs(V); 
+
+if ~isempty(pv_index)
+Theta(pv_index) = Theta(pv_index) + dx(N1:N2);
+end
+if ~isempty(pq_index)
+Theta(pq_index) = Theta(pq_index) + dx(N3:N4);
+Vm(pq_index) = Vm(pq_index) + dx(N5:N6);
+end
+
+
+V = Vm .* exp(1i * Theta);
+
+
+end

data_simulation/Power Flow Initialization/calculate_F.m

@@ -0,0 +1,8 @@
+function F = calculate_F(Ybus,Sbus,V,pv_index,pq_index)
+
+
+Delta_S = Sbus - V .* conj(Ybus * V);
+
+F = [real(Delta_S([pv_index; pq_index])); imag(Delta_S(pq_index))];
+
+end

data_simulation/Power Flow Initialization/find_equilibrium.m

@@ -0,0 +1,51 @@
+
+vq = 0;
+omega_pll =1;
+E = 1;
+Mw = omega_pll;
+
+Vref = abs(V(1));
+
+Pvsc = real(S_inj(1));
+Pref = Pvsc;
+vd = Vref;
+id = Pref/vd;
+idref = id;
+Qvsc = imag(S_inj(1));
+Qref = Qvsc;
+iq = -Qvsc/vd;
+
+
+theta_pll = atan(imag(V(1))/real((V(1))));
+
+iPcmd = id;
+iQcmq = iq;
+
+V_pcc = Vref;
+
+Vmf = V_pcc;
+E_x = E;
+E_y = 0;
+
+E_d = (E_x*cos(theta_pll));
+E_q = (-E_x*sin(theta_pll));
+
+Pext = Pref;
+Qext = Qref;
+
+
+i_x = id*cos(theta_pll) - iq*sin(theta_pll);
+i_y = id*sin(theta_pll) + iq*cos(theta_pll);
+
+P_total = E_d*id+ E_q*iq;
+Q_total = E_q*id - E_d*iq;
+
+% x0_1 = [theta_pll, Mw, id, iq, Vmf];
+% z0_1 = [vd, vq, omega_pll, Pvsc, V_pcc, Qvsc, E_d, E_q, iPcmd, iQcmq, i_x, i_y, P_total, Q_total];
+
+x0_2 = [theta_pll, Mw, id, iq, Vmf];
+z0_2 = [vd, vq, omega_pll, Pvsc, V_pcc, Qvsc, E_d, E_q, i_x, i_y, P_total, Q_total];
+
+
+
+

data_simulation/Power Flow Initialization/generate_Derivatives.m

@@ -0,0 +1,7 @@
+function [J_dS_dVm,J_dS_dTheta] = generate_Derivatives(Ybus,V)
+
+J_dS_dVm = diag(V./abs(V))*diag(conj(Ybus*V)) + diag(V)*conj(Ybus*diag(V./abs(V)));
+
+J_dS_dTheta = 1i*diag(V)*conj(diag(Ybus*V) - Ybus*diag(V));
+
+end

data_simulation/Power Flow Initialization/generate_Jacobian.m

@@ -0,0 +1,12 @@
+function J=generate_Jacobian(J_dS_dVm,J_dS_dTheta,pv_index,pq_index)
+
+
+J_11 = real(J_dS_dTheta([pv_index; pq_index], [pv_index; pq_index]));
+J_12 = real(J_dS_dVm([pv_index; pq_index], pq_index));
+J_21 = imag(J_dS_dTheta(pq_index, [pv_index; pq_index]));
+J_22 = imag(J_dS_dVm(pq_index, pq_index));
+
+J = [ J_11 J_12; J_21 J_22; ];
+
+
+end