Browse wiki

%Graphical Controls
%P_N = 4; %Number of plots (ONE OR FOUR)
%SCENARIO_NUMBERS = [ 1, 3,     4,     5];                 %DEFINE SCENARIOS TO COMPARE(between 0 and 5) IN MULTIPLOT. For single plot, first scenario is 
%V_P = [-33.75,150.5,8]; %Define map viewpoint (lat, long, zoom number). Zoom number MUST be between 6 & displayed.

L_P = [25,100]; % Legend Position ([0,0] = TOP LEFT)

for PLOTS = 1:P_N;

mpc.version = '2';

%Get the NSW network defaults baseMVA
%MachineBaseMva ?MaximumRealPowerOutput ?MinimumRealPowerOutput ?MaximumVoltageMagnitude ?MinimumVoltageMagnitude ?BusArea ?VoltageAngle ?BaseVoltage ?ShuntConductance
nswDefaults = sparql("SELECT ?MachineBaseMva ?MaximumRealPowerOutput ?MinimumRealPowerOutput ?MaximumVoltageMagnitude ?MinimumVoltageMagnitude ?BusArea ?VoltageAngle ?BaseVoltage ?ShuntConductance WHERE{
a:NSW_Electricity_Network prop:MachineBaseMva ?MachineBaseMva.
a:NSW_Electricity_Network prop:MaximumRealPowerOutput ?MaximumRealPowerOutput.
a:NSW_Electricity_Network prop:MinimumRealPowerOutput ?MinimumRealPowerOutput.
a:NSW_Electricity_Network prop:MaximumVoltageMagnitude ?MaximumVoltageMagnitude.
a:NSW_Electricity_Network prop:MinimumVoltageMagnitude ?MinimumVoltageMagnitude.
a:NSW_Electricity_Network prop:BusArea ?BusArea.
a:NSW_Electricity_Network prop:VoltageAngle ?VoltageAngle.
a:NSW_Electricity_Network prop:BaseVoltage ?BaseVoltage.
a:NSW_Electricity_Network prop:ShuntConductance ?ShuntConductance.
}");

%Parse the baseMVA
mpc.baseMVA = sscanf (char(nswDefaults(1,1)), "%f MVA");


%?bus ?id ?BusType ?realPowerDemand ?reactivPowerDemand ?shuntSusceptance ?VoltageMagnitude ?RealPowerOutput ?MaximumReactivePowerOutput ?MinimumReactivePowerOutput ?VoltageMagnitudeSetpoint ?GeoLocation ?ShortName
bus = sparql("select ?bus ?id ?BusType ?realPowerDemand ?reactivPowerDemand ?shuntSusceptance ?VoltageMagnitude ?RealPowerOutput 
  ?MaximumReactivePowerOutput ?MinimumReactivePowerOutput ?VoltageMagnitudeSetpoint ?GeoLocation ?ShortName WHERE{
?bus rdf:type cat:Electrical_Bus.
OPTIONAL{?bus prop:BusIdentifier ?id.}
OPTIONAL{?bus prop:BusType ?BusType.}
OPTIONAL{?bus prop:RealPowerDemand ?realPowerDemand.}
OPTIONAL{?bus prop:ReactivPowerDemand ?reactivPowerDemand.}
OPTIONAL{?bus prop:ShuntSusceptance ?shuntSusceptance.}
OPTIONAL{?bus prop:VoltageMagnitude ?VoltageMagnitude.}
OPTIONAL{?bus prop:RealPowerOutput ?RealPowerOutput.}
OPTIONAL{?bus prop:MaximumReactivePowerOutput ?MaximumReactivePowerOutput.}
OPTIONAL{?bus prop:MinimumReactivePowerOutput ?MinimumReactivePowerOutput.}
OPTIONAL{?bus prop:VoltageMagnitudeSetpoint ?VoltageMagnitudeSetpoint.}
OPTIONAL{?bus prop:GeoLocation ?GeoLocation.}
OPTIONAL{?bus prop:ShortName ?ShortName.}
}
ORDER BY ?bus");

[m,n] =size(bus);
B_n = n;
busCtr = 0;
genCtr = 0;

%----------------------------------------------------------------------------------------------------------------%
%Build up the mpc.bus tablefrom bus. Take defaults from the nswDefaults variable [Bus Data (mpc.bus) matrix definition ]

for i = 1:n
BUS_I(i) = cell2mat(bus(2,i)); %BUS_I  (column 1) Bus number
BUS_TYPE(i) = cell2mat(bus(3,i)); %BUS_TYPE  (colomn 2) Bus type (1 = PQ, 2 = PV, 3 = ref, 4 = isolated)
PD(i) =sscanf (char(bus(4,i)), "%d MW"); %PD  (column 3) Real Power Demand
QD(i) =sscanf (char(bus(5,i)), "%d MVA"); %QD  (column 4) Reactive Power Demand
GS(i) =  sscanf (char(nswDefaults(9,1)), "%f MW"); %GS (Column 5)  Shunt Conductance
if (length (char(bus(6,i)))<4), BS(i) = 0;  %BS (Column 6)  Shunt Susceptance
else BS(i) = sscanf (char(bus(6,i)), "%f MVA"); 
endif;
BUS_AREA(i) = cell2mat(nswDefaults(6,1)); %BUS_AREA (Column 7)  Area Number
if (length (char(bus(7,i)))<1), VM(i) = 1; 
else VM(i) = cell2mat(bus(7,i));  %VM (Column 8) Volatge Magnitude
endif;
VA(i) = cell2mat(nswDefaults(7,1)); %VA (Column 9) Volatge Angle
BASE_KV(i) = sscanf (char(nswDefaults(8,1)), "%f kV"); %BASE_KV (Column 10) Base Volatge (kV)
ZONE(i) = 1; %ZONE  (Column 11) Loss Zone
VMAX(i) = cell2mat(nswDefaults(4,1)); %VMAX (Coulmn 12)  Maximum voltage magnitude
VMIN(i) = cell2mat(nswDefaults(5,1)); %VMIN (Coulmn 13)  Minimum voltage magnitude
endfor

mpc.bus = [BUS_I', BUS_TYPE', PD',QD',GS',BS',BUS_AREA',VM',VA', BASE_KV',ZONE',VMAX',VMIN'];

%----------------------------------------------------------------------------------------------------------------%
%Build up mpc.gen table from bus. Take defaults from the nswDefaults variable

for j = 1:n

if cell2mat(bus(3,j)) ==1, GEN_BUS(j)=0; %GEN_BUS (Column 1) Bus number
else GEN_BUS(j)=cell2mat(bus(2,j)); 
endif
if (length (char(bus(8,j)))<4), PG(j) = 0;  %PG (Column 2) Real Power Output
else PG(j) = sscanf (char(bus(8,j)), "%f MW"); 
endif;
QG(j) = 0; %QG (Column 3)  Reactive Power Output  !HARD CODED
if (length (char(bus(9,j)))<4), QMAX(j) = 0;  %QMAX (Column 4) Maximum Reactive Power Output
else QMAX(j) = sscanf (char(bus(9,j)), "%f MVA"); 
endif;
if (length (char(bus(10,j)))<4), QMIN(j) = 0;  %QMIN (Column 5) Minimum Reactive Power Output
else QMIN(j) = sscanf (char(bus(10,j)), "%f MVA"); 
endif;
if (length (char(bus(11,j)))<1), VG(j) = 0;  %VG (Column 6) Voltage Magnitude
else VG(j) = cell2mat(bus(11,j)); 
endif;
MBASE(j) =  sscanf (char(nswDefaults(1,1)), "%f MVA"); %MBASE (Column 7) Total MVA base of machine !HARDCODED
GEN_STATUS(j) = 1; %GEN_STATUS (Column 8) Machine Status !HARDCODED
PMAX(j) = 9999; %PMAX (Column 9) Maximum reaal power output !HARDCODED
PMIN(j) = 9999; %PMIN (Column 10) Minumum real power output !HARDCODED
endfor

mpc.gen = [GEN_BUS', PG',QG', QMAX', QMIN', VG',MBASE',GEN_STATUS',PMAX',PMIN'];

a = mpc.gen;
[s, i] = sort (a (:, 1)); 
mpc.gen = a (i, :);
len = length (mpc.gen)-length (mpc.gen(:,1)(mpc.gen(:,1)~=0)); %Delete zero leading rows
mpc.gen(1:len,:)=[];

%----------------------------------------------------------------------------------------------------------------%
%Get the line data
%?line ?ShortName ?FromBus ?FromId ?FromGeoLocation ?ToBus ?ToId ?ToGeoLocation ?Resistance ?Reactance ?TotalLineChargingSusceptance
branch = sparql("select ?line ?ShortName ?FromBus ?FromId ?FromGeoLocation ?ToBus ?ToId ?ToGeoLocation ?Resistance ?Reactance ?TotalLineChargingSusceptance WHERE{
?line rdf:type cat:Electrical_Transmission.
?line prop:ShortName ?ShortName.
?line prop:FromBus ?FromBus.
?FromBus prop:BusIdentifier ?FromId.
?FromBus prop:GeoLocation ?FromGeoLocation.
?line prop:ToBus ?ToBus.
?ToBus prop:BusIdentifier ?ToId.
?ToBus prop:GeoLocation ?ToGeoLocation.
?line prop:Resistance ?Resistance.
?line prop:Reactance ?Reactance.
?line prop:TotalLineChargingSusceptance  ?TotalLineChargingSusceptance .
} ORDER BY ?line");

%build up the mpc.branch matrix
[m,n] =size(branch);
for i = 1:n

FBUS(i) = cell2mat(branch(4,i)); %F_Bus From Bus  (column 1)
TBUS(i) = cell2mat(branch(7,i)); %T_BUS To Bus  (column 2)
BR_R(i) = cell2mat(branch(9,i)); %BR_R Resistance  (column 3)
BR_X(i) = cell2mat(branch(10,i)); %BR_X Reactance  (column 4)
BR_B(i) = cell2mat(branch(11,i)); %BR_B Total line charging suseptance  (column 5)
RATE_A(i) = 1000; %RATE_A MVA rating A (column 6) !HARDCODED
RATE_B(i) = 1000; %RATE_B MVA rating B (column 7) !HARDCODED
RATE_C(i) = 1000; %RATE_C MVA rating C (column 8) !HARDCODED
TAP(i) = 0; %TAP Transformer off Nominal turns ratio (column 9) !HARDCODED
SHIFT(i) = 0; %SHIFT Transformer phase shift angles (column 10) !HARDCODED
BR_STATUS(i) = 1; %BR_STATUS Initial branch status (column 11) !HARDCODED
ANGMIN(i) = -360; %ANGMIN Minimum angle difference (column 12) !HARDCODED
ANGMAX(i) = 360; %ANDMAX Maximum angle difference (column 13) !HARDCODED

endfor

mpc.branch = [FBUS',TBUS',BR_R',BR_X',BR_B', RATE_A', RATE_B', RATE_C', TAP',SHIFT',BR_STATUS',ANGMIN',ANGMAX'];
mpc;

%----------------------------------------------------------------------------------------------------------------%
%Get the scenario data

%----------------------------------------------------------------------------------------------------------------%
%Evaluate the relevant case

SCENARIO = [ 1, 0.6, 1.05, 1.05, 1.05, 1.05; %SCENARIO MATRIX
1, 0.6, 1.05, 1.05, 1.05, 1.05;
1, 1, 1, 0, 1, 1;
1, 1, 1, 1, 0, 1;
1, 1, 1, 1, 1, 0 ];

SN = SCENARIO_NUMBERS(1,PLOTS)+1;

% Scenario 0 - No modification
% Scenario 1 - Base case (light load)
% Scenario 2 - Heavy load 
% Scenario 3 - Heavy loading plus outage of the lines 12-13 and 12-16 (Vales - Sydney North Vales - Munmorah)
% Scenario 4 - Heavy loading plus outage of generator #6 (Liddell)
% Scenario 5 - Heavy loading plus outage of the lines 9-19 and 9-20 (Bayswater - Regentville and Bayswater - Mt Piper)

mpc.bus(:,[3:4]) = SCENARIO(1,SN)*mpc.bus(:,[3:4]);
mpc.gen(:,[2]) = SCENARIO(2,SN)*mpc.gen(:,[2]);
mpc.branch([18,19],[11]) = SCENARIO(3,SN)*mpc.branch([18,19],[11]);
mpc.gen([2],[8]) = SCENARIO(4,SN)*mpc.gen([2],[8]);
mpc.branch([27,48],[11]) = SCENARIO(5,SN)*mpc.branch([27,48],[11]);
diary off;
results = runpf(mpc, mpoption( "OUT_ALL",0));
diary on;
%----------------------------------------------------------------------------------------------------------------%
%Apparent Power Calculation

PP_IN = diag(results.branch(:,14)*(results.branch(:,14))'); %Test Apparent power in is under limit
QQ_IN = diag(results.branch(:,15)*(results.branch(:,15))');
APPARENT_POWER_IN = sqrt (PP_IN+QQ_IN);

PP_OUT = diag(results.branch(:,16)*(results.branch(:,16))'); %Test Apparent power out is under limit
QQ_OUT = diag(results.branch(:,17)*(results.branch(:,17))');
APPARENT_POWER_OUT = sqrt (PP_OUT+QQ_OUT);

for t = 1: n
if APPARENT_POWER_IN(t,1) < 1000, APPARENT_POWER_TEST(t,1) =1; else APPARENT_POWER_TEST(t,1) = 0; endif; %Ensure BOTH power in AND power out is under limit. 
if APPARENT_POWER_OUT(t,1) < 1000, APPARENT_POWER_TEST(t,2) =1; else APPARENT_POWER_TEST(t,2) = 0; endif;
APPARENT_POWER_TEST(t,3) = APPARENT_POWER_TEST(t,2) + APPARENT_POWER_TEST(t,1);
endfor

%----------------------------------------------------------------------------------------------------------------%
%PLOTTING THE RESULTS
%----------------------------------------------------------------------------------------------------------------%


%Set Viewpoint characteristics (i.e. Lat, Long, zoom) ZOOM LEVEL MUST BE BETWEEN 6 & 8 



%V_P = [-33.75, 150.5, 8]; %Reasonable co-ordinates for zoom level 8 (Sydney level)
%V_P = [-33, 149, 7]; %Reasonable co-ordinates for zoom level 7 (North Coast level)
%V_P = [-33, 148, 6]; %Reasonable co-ordinates for zoom level 6 (NSW level)


%----------------------------------------------------------------------------------------------------------------%
%Get Image from Google Maps API

LAT = num2str(V_P(1,1));
LONG = num2str(V_P(1,2));
ZOOM = num2str(V_P(1,3));

CENTER = strcat (LAT,",",LONG);

mapbgf = urlread(strcat("http://maps.google.com/maps/api/staticmap?center=",CENTER,"&zoom=",ZOOM,"&size=640x640&sensor=false&style=feature:all",char(124),"visibility:off&style=feature:water",char(124),"saturation:-100",char(124),"lightness:-50",char(124),"visibility:on&style=feature:administrative.province",char(124),"element:geometry",char(124),"visibility:on&style=feature:landscape",char(124),"element:geometry",char(124),"lightness:100",char(124),"visibility:on"));
bgfid=fopen("static.png","w");
fwrite(bgfid,mapbgf,"char");
fclose(bgfid);

I = imread("static.png");

if P_N == 4, PLOT_GRID = 2; else PLOT_GRID = 1; endif
subplot(PLOT_GRID,PLOT_GRID,PLOTS); imshow(I)
hold on

%----------------------------------------------------------------------------------------------------------------%

%Zoom Level transform variables

T  = [ 6.05, 2.97, 1.452; %Transform matrix 
6.9, 3.5, 1.75;
45.87, 90.9, 183.15;
53.145, 108.67, 221.23 ];

Z_n = V_P(1,3) - 5; %Zoom Number (1,2 or 3, where 1 = "zoom level 6", 2 = "zoom level 7", 3 = "zoom level 8")


Y = -V_P (1,1) - T(1,Z_n); %translation in y
X = V_P(1,2) - T(2,Z_n); %translation in x
XD = T(3,Z_n); %dilation factor x
YD = T(4,Z_n); %dilation factor y


%----------------------------------------------------------------------------------------------------------------%

%Color scheme


%Ok lines- 4284d3 – 66,132,211, width 1 
LINES = [11,97,164]/256;
%Err lines - FFB540 – 255,181,64, width 2
ERRLINES = [255,191,0]/256;
%Out lines – 000000 – 0,0,0 – width 1, could we make these dashed?
OUTLINES = [255,49,0]/256;
%Ok bus - 04346C – 4, 52,108 – Radius 1 - Dot
BUS_COL = [3,62,107]/256;
%Err bus - A66500 – 166,101,0 – Radius 2 – Concentric circles
ERR_BUS_COL = [166,124,0]/256;

%----------------------------------------------------------------------------------------------------------------%
subplot(PLOT_GRID,PLOT_GRID,PLOTS);

%Plot Branches
for i = 1:n
c = [sscanf(char(branch(5,i)), '%f,'),sscanf(char(branch(8,i)), '%f,')];
if ((APPARENT_POWER_IN(i,1)) == 0),
plot (((c(2,:)-X)*XD),((Y+c(1,:))*-YD),"color",OUTLINES,"linewidth",2,"linestyle","--")
elseif (abs (APPARENT_POWER_TEST(i,3)) < 2), 
plot (((c(2,:)-X)*XD),((Y+c(1,:))*-YD),"color",ERRLINES,"linewidth",2)
else plot (((c(2,:)-X)*XD),((Y+c(1,:))*-YD),"color",LINES,"linewidth",1)
endif
endfor


%Plot Buses
for i = 1:B_n
GL = sscanf(char(bus(12,i)), '%f,');
if results.bus(i,8) > 1.1
for k = 1:3 plot (((GL(2,1)-X)*XD),((Y+GL(1,1))*-YD), "o","color",ERR_BUS_COL, "markersize",(2*k)) endfor
elseif results.bus(i,8) < 0.9 
for k = 1:3 plot (((GL(2,1)-X)*XD),((Y+GL(1,1))*-YD), "o","color",ERR_BUS_COL,"markersize",(2*k)) endfor
else for k = 1:3 plot (((GL(2,1)-X)*XD),((Y+GL(1,1))*-YD), "o","color",BUS_COL,"markersize",k) endfor
endif
endfor

%----------------------------------------------------------------------------------------------------------------%

%Plot Legend



L_L = 3; % Line length in legend
if P_N == 4, L_S = 20; else L_S = 10; endif; % Legend Spacing  
if P_N == 4, L_TS = 30; else L_TS = 15; endif; % Legend Title spacing

%Legend images
for k = 1:3 plot (L_P(1,1), L_P(1,2)+L_TS, "o","color",ERR_BUS_COL,"markersize",(2*k)) endfor
for k = 1:3 plot (L_P(1,1), L_P(1,2)+L_TS+L_S, "o","color",BUS_COL,"markersize",k) endfor
plot ([(L_P(1,1)-L_L),(L_P(1,1)+L_L)],[L_P(1,2)+L_TS+2*L_S,L_P(1,2)+L_TS+2*L_S],"color",OUTLINES,"linewidth",2)
plot ([(L_P(1,1)-L_L),(L_P(1,1)+L_L)],[L_P(1,2)+L_TS+3*L_S,L_P(1,2)+L_TS+3*L_S],"color",ERRLINES,"linewidth",2)
plot ([(L_P(1,1)-L_L),(L_P(1,1)+L_L)],[L_P(1,2)+L_TS+4*L_S,L_P(1,2)+L_TS+4*L_S],"color",LINES,"linewidth",1)

%Legend text
text (L_P(1,1), L_P(1,2), "LEGEND")
text (L_P(1,1)+10, L_P(1,2)+L_TS, "Bus Voltage Outside Limits")
text (L_P(1,1)+10, L_P(1,2)+L_TS+L_S, "Bus OK")
text (L_P(1,1)+10, L_P(1,2)+L_TS+2*L_S, "Line Outage")
text (L_P(1,1)+10, L_P(1,2)+L_TS+3*L_S, "Line Overload")
text (L_P(1,1)+10, L_P(1,2)+L_TS+4*L_S, "Line Ok")

%TITLE

SCEN_TITLE = strcat ("Scenario",{" "},num2str(SN-1));
text (300,0, SCEN_TITLE)

%----------------------------------------------------------------------------------------------------------------%
endfor

hold off