-
Notifications
You must be signed in to change notification settings - Fork 7
/
DistFlowObjectUnstable.m
111 lines (93 loc) · 2.67 KB
/
DistFlowObjectUnstable.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
clc
clear
close all
FromNode=[1 2 3 4 5];
ToNode=[2 3 4 5 6];
Branch=1:5;
%% Volt Var Volt Watt Curves
linear=[0.95,0.98,1.02,1.05,1.02,1.05];
% linear=linear';
no_deadband=[0.95,1.00,1.00,1.05,1.0,1.05];
curved=[0.95,0.98,1.02,1.05,1.05,1.1];
gen_control=[0,0,1,0,0,1];
% gen_control=gen_control*0;
Smaxa=[0,40,20,10,30,50];
% Pmax=MaxCollection(Smaxa);
%%
F=sparse(Branch,FromNode,ones(1,5),5,6);
T=sparse(Branch,ToNode,ones(1,5),5,6);
M=F-T;
r=0.01;
x=0.01;
r=r*[1 1.2 1.2 1.3 0.95];
x=x*[1 1.2 1.2 1.3 0.95];
NumberOfBranch=5;
NumberOfNodes=6;
Rline=r.*speye(NumberOfBranch,NumberOfBranch);
Xline=x.*speye(NumberOfBranch,NumberOfBranch);
TFT=T*F';
I=sparse(eye(size(TFT)));
originalM=M;
M=M(:,2:end);
R=2*(M\Rline)*((I-TFT)\T);
X=2*(M\Xline)*((I-TFT)\T);
slack_voltage=1.01;
% load values
% pc=[0,0.3,0.9,0.2,0,0.3];
% qc=[0,0.8,0.2,0.3,0.1,0.2];
pc=[0,0.3,0.0,0.2,0.9,0.3];
qc=[0,0.0,0.2,0,0.1,0.2];
%% Inverter Specifications
% Number_of_Inverters=4;
% stable case
% InverterSmax=[20,50,40];
InverterNode=[2,3];
% InverterNode=[3,5,6];
InverterSmax=[40,40];
% droop_curves={'no_deadband','no_deadband'};
droop_curves={'no_deadband','no_deadband','curved'};
droop=containers.Map;
droop('curved')=curved;
droop('linear')=linear;
droop('no_deadband')=no_deadband;
Number_of_Inverters=length(InverterNode);
InverterArray(1:Number_of_Inverters)=InverterDistFlow;
for i = 1:Number_of_Inverters
InverterArray(i).Smax=InverterSmax(i);
InverterArray(i)=MaxCollection(InverterArray(i));
InverterArray(i).Node=InverterNode(i);
InverterArray(i).setpoints=droop(droop_curves{i});
end
%% Unstable Case
%%
iteration=300;
VoltageArray=[];
Pg=[];
Qg=[];
for itr=1:iteration
if (itr==1)
% Solving the first iteration without inverter
V=[slack_voltage^2 ;slack_voltage^2+R*pc'+X*qc'];
VoltageArray=V';
Pg=zeros(1,NumberOfNodes);
Qg=zeros(1,NumberOfNodes);
else
pg=zeros(1,NumberOfNodes);
qg=zeros(1,NumberOfNodes);
% Start Iteration for All inverters
for i = 1:Number_of_Inverters
pg(1,InverterArray(i).Node)=PCalc(InverterArray(i),V(InverterArray(i).Node));
qg(1,InverterArray(i).Node)=QCalc(InverterArray(i),V(InverterArray(i).Node));
end
Pg=[Pg; pg];
Qg=[Qg; qg];
V=[slack_voltage^2 ; slack_voltage^2+R*(pc'-0.01*pg')+X*(qc'-0.01*qg')];
VoltageArray=[VoltageArray;V'];
diffv=VoltageArray(itr,:)-VoltageArray(itr-1,:);
diffv=abs(diffv);
if max(diffv) < 0.0001
break;
end
end
end
% disp(sqrt(VoltageArray))