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 define_constants;
0047
0048 if nargin < 2
0049 mpopt = mpoption;
0050 end
0051
0052
0053 ignore_angle_lim = mpopt.opf.ignore_angle_lim;
0054 verbose = mpopt.verbose;
0055
0056 if ~have_feature('yalmip')
0057 error('insolvablepfsos_limitQ: The software package YALMIP is required to run insolvablepfsos_limitQ. See https://yalmip.github.io');
0058 end
0059
0060
0061 sdpopts = yalmip_options([], mpopt);
0062
0063
0064
0065 mpc = loadcase(mpc);
0066 mpc = ext2int(mpc);
0067
0068 if toggle_dcline(mpc, 'status')
0069 error('insolvablepfsos_limitQ: DC lines are not implemented in insolvablepf');
0070 end
0071
0072 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360))
0073 warning('insolvablepfsos_limitQ: Angle difference constraints are not implemented. Ignoring angle difference constraints.');
0074 end
0075
0076 nbus = size(mpc.bus,1);
0077 ngen = size(mpc.gen,1);
0078
0079
0080
0081
0082 for i=1:ngen
0083 busidx = find(mpc.bus(:,BUS_I) == mpc.gen(i,GEN_BUS));
0084 if isempty(busidx) || ~(mpc.bus(busidx,BUS_TYPE) == PV || mpc.bus(busidx,BUS_TYPE) == REF)
0085 mpc.bus(busidx,BUS_TYPE) = PV;
0086 if verbose >= 1
0087 warning('insolvablepfsos_limitQ: Bus %s has generator(s) but is listed as a PQ bus. Changing to a PV bus.',int2str(busidx));
0088 end
0089 end
0090 end
0091
0092
0093
0094 for i=1:nbus
0095 if mpc.bus(i,BUS_TYPE) == PV
0096 genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I), 1);
0097 if isempty(genidx)
0098 mpc.bus(i,BUS_TYPE) = PQ;
0099 if verbose >= 1
0100 warning('insolvablepfsos_limitQ: PV bus %i has no associated generator! Changing these buses to PQ.',i);
0101 end
0102 end
0103 end
0104 end
0105
0106 pq = find(mpc.bus(:,2) == PQ);
0107 pv = find(mpc.bus(:,2) == PV);
0108 ref = find(mpc.bus(:,2) == REF);
0109
0110 refpv = sort([ref; pv]);
0111 pvpq = sort([pv; pq]);
0112
0113 Y = makeYbus(mpc);
0114 G = real(Y);
0115 B = imag(Y);
0116
0117
0118 S = makeSbus(mpc.baseMVA, mpc.bus, mpc.gen);
0119 P = real(S);
0120 Q = imag(S);
0121
0122 Vmag = zeros(nbus,1);
0123 Vmag(mpc.gen(:,GEN_BUS)) = mpc.gen(:,VG);
0124
0125
0126 Qmax = 1e3*ones(length(refpv),1);
0127 for i=1:length(refpv)
0128 genidx = find(mpc.gen(:,1) == refpv(i));
0129 Qmax(i) = (sum(mpc.gen(genidx,QMAX)) - mpc.bus(mpc.gen(genidx(1),1),QD)) / mpc.baseMVA;
0130
0131 end
0132
0133
0134
0135
0136
0137
0138 if verbose >= 2
0139 fprintf('Creating variables.\n');
0140 end
0141
0142 Vd = sdpvar(nbus,1);
0143 Vq = sdpvar(nbus,1);
0144
0145
0146
0147 Vminus2 = sdpvar(ngen,1);
0148 x = sdpvar(ngen,1);
0149
0150
0151 V = [1; Vd; Vq; Vminus2; x];
0152
0153
0154
0155 if verbose >= 2
0156 fprintf('Creating polynomials.\n');
0157 end
0158
0159 deg = 0;
0160 sosdeg = 0;
0161
0162
0163 for i=1:nbus-1
0164 [l,lc] = polynomial(V,deg,0);
0165 lambda(i) = l;
0166 clambda{i} = lc;
0167 end
0168
0169
0170 for i=1:length(pq)
0171 [g,gc] = polynomial(V,deg,0);
0172 gamma(i) = g;
0173 cgamma{i} = gc;
0174 end
0175
0176
0177 [delta, cdelta] = polynomial(V,deg,0);
0178
0179
0180 for i=1:length(refpv)
0181 [z,cz] = polynomial(V,deg,0);
0182 zeta{i} = z;
0183 zeta{i}(1) = z;
0184 czeta{i}{1} = cz;
0185
0186 [z,cz] = polynomial(V,deg,0);
0187 zeta{i}(2) = z;
0188 czeta{i}{2} = cz;
0189
0190 [z,cz] = polynomial(V,deg,0);
0191 zeta{i}(3) = z;
0192 czeta{i}{3} = cz;
0193
0194
0195
0196
0197
0198 [s1,c1] = polynomial(V,sosdeg,0);
0199 [s2,c2] = polynomial(V,sosdeg,0);
0200
0201
0202
0203 sigma{i} = s1;
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 if have_feature('matlab', 'vnum') >= 8.006 && have_feature('cplex') && ...
0366 have_feature('cplex', 'vnum') <= 12.006003
0367 warning('OFF', 'MATLAB:lang:badlyScopedReturnValue');
0368 end
0369
0370 sol = solvesos(constraints,[],sosopts,pvec);
0371
0372 if ~have_feature('octave') || have_feature('octave', 'vnum') >= 4.001
0373
0374 warning(S);
0375 end
0376
0377 if verbose >= 2
0378 fprintf('Solver exit message: %s\n',sol.info);
0379 end
0380
0381 if sol.problem
0382
0383 insolvable = 0;
0384 elseif sol.problem == 0
0385
0386 insolvable = 1;
0387 end