-
Notifications
You must be signed in to change notification settings - Fork 11
/
runPF2
143 lines (109 loc) · 5.86 KB
/
runPF2
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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
function [V, I_1, I_2] = runPF2(terminal2Node, DSSObj)
% FUNCTION runPF: perform steady state power flow calculation and return
% complex bus voltage and current
% INPUT:
% mappingTerminal2Bus: bidirectional mapping between nodes and buses
% (including phase identifier)
% DSSObj: a DSS Object
% OUTPUT:
% V: a vector of size NODES containing complex voltage at each node
% I_1: a vector of size NODES containing sending-end complex current at each node
% I_2: a vector of size NODES containing receiving-end complex current at each node
% Define the text interface
DSSText = DSSObj.Text;
%Reference the circuit for the interface
DSSCircuit = DSSObj.ActiveCircuit;
DSSSolution= DSSCircuit.Solution;
nNodes = length(terminal2Node.keys);
nodeList = terminal2Node.keys;
phaseIndex = ['a', 'b', 'c'];
LoadBusNamesArray = [];
LoadBusRealLoadArray = [];
LoadBusReactiveLoadArray = [];
lelem = DSSCircuit.Loads.First;
while lelem>0
LoadBusNamesArray = [LoadBusNamesArray,cellstr(DSSCircuit.Loads.Name)];
LoadBusRealLoadArray = [LoadBusRealLoadArray;num2cell(DSSCircuit.Loads.kW)];
LoadBusReactiveLoadArray = [LoadBusReactiveLoadArray;num2cell(DSSCircuit.Loads.kvar)];
% PL = DSSCircuit.Loads.kW+normrnd(0,.01*DSSCircuit.Loads.kW);
% QL = DSSCircuit.Loads.kvar+normrnd(0,.01*DSSCircuit.Loads.kvar);
PL = normrnd(7,1);
QL = normrnd(5,.5);
DSSText.Command = ['Edit Load.' DSSCircuit.Loads.Name ' kW=' num2str(PL) ' kvar=' num2str(QL)];
lelem = DSSCircuit.Loads.Next;
end
DSSText.Command = 'vsource.source.enabled = yes';
%DSSText.Command = 'Clear';
%DSSText.Command='Solve mode=snap'
% perform power flow calculation
DSSSolution.Solve;
for i=1:length(LoadBusNamesArray)
DSSText.Command = ['Edit Load.' char(LoadBusNamesArray(i)) ' kW=' num2str(cell2mat(LoadBusRealLoadArray(i))) ' kvar=' num2str(cell2mat(LoadBusReactiveLoadArray(i)))];
end
% check if a solution is found
if DSSSolution.Converged
% get bus voltages
vckt = DSSCircuit.YNodeVarray;
vckt_len = length(vckt);
V = vckt(1:2:vckt_len)+1i*vckt(2:2:vckt_len);
% sender-end currents
I_1 = zeros(1,nNodes);
% receiving-end currents
I_2 = zeros(1,nNodes);
lelem = DSSCircuit.Lines.First;
while lelem>0
DeviceIndex = DSSCircuit.SetActiveElement(['Line.' DSSCircuit.Lines.Name]);
% get buses that are connected by a given line
blist = DSSCircuit.CktElements(int32(DeviceIndex)).BusNames;
%display(blist)
% extract the sender node name
senderBus = strtok(blist{1},'.');
% extract the receiver node name
receiverBus = strtok(blist{2},'.');
% extract currents passing through the line (real and img components)
CArray = DSSCircuit.CktElements(int32(DeviceIndex)).Currents;
% this returns real and img parts of current flow for each phase in A (once from sender to receiver then from receiver to sender)
OArray = DSSCircuit.CktElements(int32(DeviceIndex)).NodeOrder;
for p=1:length(CArray)/4
% get complex current at the sending-end
I_1(terminal2Node([senderBus '.' num2str(OArray(p))])) = I_1(terminal2Node([senderBus '.' num2str(OArray(p))])) + CArray(2*p-1) + CArray(2*p)*1j;
% display(['inserting ' num2str(CArray(2*p-1)) '+j' num2str(CArray(2*p)) ' into I1(' int2str(terminal2Node([senderBus '.' num2str(OArray(p))])) '): bus ' senderBus])
end
for p=1+length(CArray)/4:length(CArray)/2
% write complex current at the receiving-end
I_2(terminal2Node([receiverBus '.' num2str(OArray(p))])) = I_2(terminal2Node([receiverBus '.' num2str(OArray(p))])) + CArray(2*p-1) + CArray(2*p)*1j;
% display(['inserting ' num2str(CArray(2*p-1)) '+j' num2str(CArray(2*p)) ' into I2(' int2str(terminal2Node([receiverBus '.' num2str(OArray(p))])) '): bus ' receiverBus])
end
lelem = DSSCircuit.Lines.Next;
end
telem = DSSCircuit.Transformers.First;
while telem>0
DeviceIndex = DSSCircuit.SetActiveElement(['Transformer.' DSSCircuit.Transformers.Name]);
% get buses that are connected by a given line
blist = DSSCircuit.CktElements(int32(DeviceIndex)).BusNames;
%display(blist)
% extract the sender node name
senderBus = strtok(blist{1},'.');
% extract the receiver node name
receiverBus = strtok(blist{2},'.');
% extract currents passing through the line (real and img components)
CArray = DSSCircuit.CktElements(int32(DeviceIndex)).Currents;
%display(CArray)
% this returns real and img parts of current flow for each phase in A (once from sender to receiver then from receiver to sender)
OArray = DSSCircuit.CktElements(int32(DeviceIndex)).NodeOrder;
for p=1:length(CArray)/4-1
% get complex current at the sending-end
I_1(terminal2Node([senderBus '.' num2str(OArray(p))])) = I_1(terminal2Node([senderBus '.' num2str(OArray(p))])) + CArray(2*p-1) + CArray(2*p)*1j;
% display(['inserting ' num2str(CArray(2*p-1)) '+j' num2str(CArray(2*p)) ' into I1(' int2str(terminal2Node([senderBus '.' num2str(OArray(p))])) '): bus ' senderBus])
end
for p=1+length(CArray)/4:length(CArray)/2-1
% write complex current at the receiving-end
I_2(terminal2Node([receiverBus '.' num2str(OArray(p))])) = I_2(terminal2Node([receiverBus '.' num2str(OArray(p))])) + CArray(2*p-1) + CArray(2*p)*1j;
% display(['inserting ' num2str(CArray(2*p-1)) '+j' num2str(CArray(2*p)) ' into I2(' int2str(terminal2Node([receiverBus '.' num2str(OArray(p))])) '): bus ' receiverBus])
end
telem = DSSCircuit.Transformers.Next;
end
else
display('did not converge')
error('Fatal error!\n OpenDSS was unable to solve PF!');
end