-
Notifications
You must be signed in to change notification settings - Fork 0
/
myHUnbalanced.m
124 lines (105 loc) · 2.63 KB
/
myHUnbalanced.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
112
113
114
115
116
117
118
119
120
121
122
123
function [Cs,bs,Cr,br,H,hT0source] = myHUnbalanced(mh,lh,sh,KT,HTssourse,hT0,HTrload,Cp,mpc)
% function [Cs,bs] = myHUnbalanced(mpc,KT,Cp,Tssourse)
[FBUS, TBUS, LENGTH, DIAMETER, FBRANCH] = myH_idx_branch;
[LOAD, SOURSE, NONE, BUS_I, BUS_TYPE, PHID, TO, TS, TR] = myH_idx_bus;
bus = mpc.bus;
branch = mpc.branch;
nb = length(bus);
nl = length(branch);
[load,sourse] = myH_bustypes(bus);
%计算hk
H = exp(-KT*branch(:,LENGTH)./(branch(:,FBRANCH)*Cp));
%% 计算Cs,bs
f = branch(:, FBUS); %% list of "from" buses
t = branch(:, TBUS); %% list of "to" buses
Cs = sparse(t,f,-H.*branch(:,FBRANCH),nb,nb);
Cs(sourse,:) = [];
bs = -sum(Cs(:,sourse),2).*HTssourse;
Cs(:,sourse) = [];
Csdiag = sum(sparse(t,f,branch(:,FBRANCH),nb,nb),2);
Csdiag(sourse,:) = [];
Cs = Cs + diag(Csdiag);
%% 计算负荷节点Cr,br
Cr=zeros(lh,lh);
for x=1:lh
% disp('节点号');
% disp(x);
a=0;
b=1;
t=1;
sss=zeros(1,mh);
rrr=zeros(1,mh);
for i=1:mh
if branch(i,FBUS)==x
a=a+1;
rrr(b)=i;%把流入i的支路序号存入rrr
ppp(b)=branch(i,2);
% disp('a');
% disp(a);
% disp(rrr);
t=b;
b=b+1;
end
end
c=1;
w=1;
for i=1:mh
if branch(i,2)==x
sss(c)=i;%把流出i的支路序号存入sss
% disp(sss);
w=c;
c=c+1;
end
end
%%类型一:只有一个负荷流回节点x
if a==0
br(x,1)=hT0(x);
Cr(x,x)=1;
end
%%类型二:负荷+多条管道流回该节点x
if a>0
bt=0;
cr=0;
for i=1:t
bt=bt-branch(rrr(i),FBRANCH);
Cr(x,ppp(i))=-branch(rrr(i),FBRANCH)*H(rrr(i));
end
for i=1:w
cr=cr+branch(sss(i),FBRANCH);
end
br(x,1)=(cr+bt)*hT0(x);
Cr(x,x)= cr;
end
end
% disp('Cr:');
% disp(Cr);
% disp('br:');
% disp(br);
%%计算源处回热温度
for x=1:sh
a=0;
b=1;
t=1;
for i=1:mh
if branch(i,FBUS)==(lh+x)
a=a+1;
mmm(b)=i;%把与源相连的支路号存入mmm
nnn(b)=branch(i,TBUS);%把与源相连的节点号存入mmm
t=b;
b=b+1;
end
end
if a==1
hT0source(x)=H(mmm(1))*HTrload(nnn(1));
end
if a>1
fenzi=0;
fenmu=0;
for i=1:t
fenzi=fenzi+H(mmm(i))*branch(mmm(i),FBRANCH)*HTrload(nnn(i));
fenmu=fenmu+branch(mmm(i),FBRANCH);
end
hT0source(x)=fenzi/fenmu;
end
end
end