0001 function insolvable = insolvablepfsos_limitQ(mpc,mpopt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067 define_constants;
0068
0069 if nargin < 2
0070 mpopt = mpoption;
0071 end
0072
0073
0074 ignore_angle_lim = mpopt.opf.ignore_angle_lim;
0075 verbose = mpopt.verbose;
0076
0077 if ~have_fcn('yalmip')
0078 error('insolvablepfsos_limitQ: The software package YALMIP is required to run insolvablepfsos_limitQ. See http://users.isy.liu.se/johanl/yalmip/');
0079 end
0080
0081
0082 sdpopts = yalmip_options([], mpopt);
0083
0084
0085
0086 mpc = loadcase(mpc);
0087 mpc = ext2int(mpc);
0088
0089 if toggle_dcline(mpc, 'status')
0090 error('insolvablepfsos_limitQ: DC lines are not implemented in insolvablepf');
0091 end
0092
0093 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360))
0094 warning('insolvablepfsos_limitQ: Angle difference constraints are not implemented. Ignoring angle difference constraints.');
0095 end
0096
0097 nbus = size(mpc.bus,1);
0098 ngen = size(mpc.gen,1);
0099
0100
0101
0102
0103 for i=1:ngen
0104 busidx = find(mpc.bus(:,BUS_I) == mpc.gen(i,GEN_BUS));
0105 if isempty(busidx) || ~(mpc.bus(busidx,BUS_TYPE) == PV || mpc.bus(busidx,BUS_TYPE) == REF)
0106 mpc.bus(busidx,BUS_TYPE) = PV;
0107 if verbose >= 1
0108 warning('insolvablepfsos_limitQ: Bus %s has generator(s) but is listed as a PQ bus. Changing to a PV bus.',int2str(busidx));
0109 end
0110 end
0111 end
0112
0113
0114
0115 for i=1:nbus
0116 if mpc.bus(i,BUS_TYPE) == PV
0117 genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I), 1);
0118 if isempty(genidx)
0119 mpc.bus(i,BUS_TYPE) = PQ;
0120 if verbose >= 1
0121 warning('insolvablepfsos_limitQ: PV bus %i has no associated generator! Changing these buses to PQ.',i);
0122 end
0123 end
0124 end
0125 end
0126
0127 pq = find(mpc.bus(:,2) == PQ);
0128 pv = find(mpc.bus(:,2) == PV);
0129 ref = find(mpc.bus(:,2) == REF);
0130
0131 refpv = sort([ref; pv]);
0132 pvpq = sort([pv; pq]);
0133
0134 Y = makeYbus(mpc);
0135 G = real(Y);
0136 B = imag(Y);
0137
0138
0139 S = makeSbus(mpc.baseMVA, mpc.bus, mpc.gen);
0140 P = real(S);
0141 Q = imag(S);
0142
0143 Vmag = zeros(nbus,1);
0144 Vmag(mpc.gen(:,GEN_BUS)) = mpc.gen(:,VG);
0145
0146
0147 Qmax = 1e3*ones(length(refpv),1);
0148 for i=1:length(refpv)
0149 genidx = find(mpc.gen(:,1) == refpv(i));
0150 Qmax(i) = (sum(mpc.gen(genidx,QMAX)) - mpc.bus(mpc.gen(genidx(1),1),QD)) / mpc.baseMVA;
0151
0152 end
0153
0154
0155
0156
0157
0158
0159 if verbose >= 2
0160 fprintf('Creating variables.\n');
0161 end
0162
0163 Vd = sdpvar(nbus,1);
0164 Vq = sdpvar(nbus,1);
0165
0166
0167
0168 Vminus2 = sdpvar(ngen,1);
0169 x = sdpvar(ngen,1);
0170
0171
0172 V = [1; Vd; Vq; Vminus2; x];
0173
0174
0175
0176 if verbose >= 2
0177 fprintf('Creating polynomials.\n');
0178 end
0179
0180 deg = 0;
0181 sosdeg = 0;
0182
0183
0184 for i=1:nbus-1
0185 [l,lc] = polynomial(V,deg,0);
0186 lambda(i) = l;
0187 clambda{i} = lc;
0188 end
0189
0190
0191 for i=1:length(pq)
0192 [g,gc] = polynomial(V,deg,0);
0193 gamma(i) = g;
0194 cgamma{i} = gc;
0195 end
0196
0197
0198 [delta, cdelta] = polynomial(V,deg,0);
0199
0200
0201 for i=1:length(refpv)
0202 [z,cz] = polynomial(V,deg,0);
0203 zeta{i}(1) = z;
0204 czeta{i}{1} = cz;
0205
0206 [z,cz] = polynomial(V,deg,0);
0207 zeta{i}(2) = z;
0208 czeta{i}{2} = cz;
0209
0210 [z,cz] = polynomial(V,deg,0);
0211 zeta{i}(3) = z;
0212 czeta{i}{3} = cz;
0213
0214
0215
0216
0217
0218 [s1,c1] = polynomial(V,sosdeg,0);
0219 [s2,c2] = polynomial(V,sosdeg,0);
0220
0221
0222
0223 sigma{i}(1) = s1;
0224 csigma{i}{1} = c1;
0225 sigma{i}(2) = s2;
0226 csigma{i}{2} = c2;
0227
0228
0229
0230
0231
0232
0233 end
0234
0235
0236
0237 if verbose >= 2
0238 fprintf('Creating power flow equations.\n');
0239 end
0240
0241
0242 for k=1:nbus-1
0243 i = pvpq(k);
0244
0245
0246
0247
0248 y = find(G(:,i) | B(:,i));
0249 for m=1:length(y)
0250 if exist('fp','var') && length(fp) == k
0251 fp(k) = fp(k) + Vd(i) * (G(y(m),i)*Vd(y(m)) - B(y(m),i)*Vq(y(m))) + Vq(i) * (B(y(m),i)*Vd(y(m)) + G(y(m),i)*Vq(y(m)));
0252 else
0253 fp(k) = Vd(i) * (G(y(m),i)*Vd(y(m)) - B(y(m),i)*Vq(y(m))) + Vq(i) * (B(y(m),i)*Vd(y(m)) + G(y(m),i)*Vq(y(m)));
0254 end
0255 end
0256
0257 end
0258
0259
0260 for k=1:length(pq)
0261 i = pq(k);
0262
0263
0264
0265
0266 y = find(G(:,i) | B(:,i));
0267 for m=1:length(y)
0268 if exist('fq','var') && length(fq) == k
0269 fq(k) = fq(k) + Vd(i) * (-B(y(m),i)*Vd(y(m)) - G(y(m),i)*Vq(y(m))) + Vq(i) * (G(y(m),i)*Vd(y(m)) - B(y(m),i)*Vq(y(m)));
0270 else
0271 fq(k) = Vd(i) * (-B(y(m),i)*Vd(y(m)) - G(y(m),i)*Vq(y(m))) + Vq(i) * (G(y(m),i)*Vd(y(m)) - B(y(m),i)*Vq(y(m)));
0272 end
0273 end
0274 end
0275
0276
0277 for k=1:length(refpv)
0278 i = refpv(k);
0279 fv(k) = Vd(i)^2 + Vq(i)^2;
0280 end
0281
0282
0283 for k=1:length(refpv)
0284 i = refpv(k);
0285
0286
0287
0288
0289 y = find(G(:,i) | B(:,i));
0290 for m=1:length(y)
0291 if exist('fq_gen','var') && length(fq_gen) == k
0292 fq_gen(k) = fq_gen(k) + Vd(i) * (-B(y(m),i)*Vd(y(m)) - G(y(m),i)*Vq(y(m))) + Vq(i) * (G(y(m),i)*Vd(y(m)) - B(y(m),i)*Vq(y(m)));
0293 else
0294 fq_gen(k) = Vd(i) * (-B(y(m),i)*Vd(y(m)) - G(y(m),i)*Vq(y(m))) + Vq(i) * (G(y(m),i)*Vd(y(m)) - B(y(m),i)*Vq(y(m)));
0295 end
0296 end
0297 end
0298
0299
0300
0301
0302 if verbose >= 2
0303 fprintf('Forming the test polynomial.\n');
0304 end
0305
0306 sospoly = 1;
0307
0308
0309 for i=1:nbus-1
0310 sospoly = sospoly + lambda(i) * (fp(i) - P(pvpq(i)));
0311 end
0312
0313
0314 for i=1:length(pq)
0315 sospoly = sospoly + gamma(i) * (fq(i) - Q(pq(i)));
0316 end
0317
0318
0319 sospoly = sospoly + delta * Vq(ref);
0320
0321
0322 for k=1:length(refpv)
0323 i = refpv(k);
0324
0325
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337
0338 sospoly = sospoly + zeta{k}(1) * (Vmag(i)^2 - Vminus2(k) - fv(k));
0339 sospoly = sospoly + zeta{k}(2) * (Vminus2(k) * x(k));
0340 sospoly = sospoly + zeta{k}(3) * (Qmax(k) - x(k) - fq_gen(k));
0341
0342 sospoly = sospoly + sigma{k}(1) * Vminus2(k);
0343 sospoly = sospoly + sigma{k}(2) * x(k);
0344 end
0345
0346 sospoly = -sospoly;
0347
0348
0349
0350 if verbose >= 2
0351 fprintf('Solving the sum of squares problem.\n');
0352 end
0353
0354 sosopts = sdpsettings(sdpopts,'sos.newton',1,'sos.congruence',1);
0355
0356 constraints = sos(sospoly);
0357
0358 pvec = cdelta(:).';
0359
0360 for i=1:length(lambda)
0361 pvec = [pvec clambda{i}(:).'];
0362 end
0363
0364 for i=1:length(gamma)
0365 pvec = [pvec cgamma{i}(:).'];
0366 end
0367
0368 for i=1:length(refpv)
0369 for k=1:length(czeta{i})
0370 pvec = [pvec czeta{i}{k}(:).'];
0371 end
0372
0373 for k=1:length(csigma{i})
0374 pvec = [pvec csigma{i}{k}(:).'];
0375 end
0376
0377 for k=1:length(csigma{i})
0378 constraints = [constraints; sos(sigma{i}(k))];
0379 end
0380 end
0381
0382
0383 S = warning;
0384
0385 sol = solvesos(constraints,[],sosopts,pvec);
0386
0387 warning(S);
0388
0389 if verbose >= 2
0390 fprintf('Solver exit message: %s\n',sol.info);
0391 end
0392
0393 if sol.problem
0394
0395 insolvable = 0;
0396 elseif sol.problem == 0
0397
0398 insolvable = 1;
0399 end