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