0001 function insolvable = insolvablepfsos(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 define_constants;
0042
0043 if nargin < 2
0044 mpopt = mpoption;
0045 end
0046
0047
0048 ignore_angle_lim = mpopt.opf.ignore_angle_lim;
0049 verbose = mpopt.verbose;
0050 enforce_Qlimits = mpopt.pf.enforce_q_lims;
0051
0052 if enforce_Qlimits > 0
0053 enforce_Qlimits = 1;
0054 end
0055
0056 if ~have_fcn('yalmip')
0057 error('insolvablepfsos: The software package YALMIP is required to run insolvablepfsos. See http://users.isy.liu.se/johanl/yalmip/');
0058 end
0059
0060
0061 sdpopts = yalmip_options([], mpopt);
0062
0063
0064
0065
0066
0067
0068 if enforce_Qlimits
0069 if verbose > 0
0070 fprintf('Generator reactive power limits are enforced. Using function insolvablepfsos_limitQ.\n');
0071 end
0072 insolvable = insolvablepfsos_limitQ(mpc,mpopt);
0073 return;
0074 end
0075
0076
0077
0078 mpc = loadcase(mpc);
0079 mpc = ext2int(mpc);
0080
0081 if toggle_dcline(mpc, 'status')
0082 error('insolvablepfsos: DC lines are not implemented in insolvablepf');
0083 end
0084
0085 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360))
0086 warning('insolvablepfsos: Angle difference constraints are not implemented. Ignoring angle difference constraints.');
0087 end
0088
0089 nbus = size(mpc.bus,1);
0090 ngen = size(mpc.gen,1);
0091
0092
0093
0094
0095 for i=1:ngen
0096 busidx = find(mpc.bus(:,BUS_I) == mpc.gen(i,GEN_BUS));
0097 if isempty(busidx) || ~(mpc.bus(busidx,BUS_TYPE) == PV || mpc.bus(busidx,BUS_TYPE) == REF)
0098 mpc.bus(busidx,BUS_TYPE) = PV;
0099 if verbose >= 1
0100 warning('insolvablepfsos: Bus %s has generator(s) but is listed as a PQ bus. Changing to a PV bus.',int2str(busidx));
0101 end
0102 end
0103 end
0104
0105
0106
0107 for i=1:nbus
0108 if mpc.bus(i,BUS_TYPE) == PV
0109 genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I), 1);
0110 if isempty(genidx)
0111 mpc.bus(i,BUS_TYPE) = PQ;
0112 if verbose >= 1
0113 warning('insolvablepfsos: PV bus %i has no associated generator! Changing these buses to PQ.',i);
0114 end
0115 end
0116 end
0117 end
0118
0119 pq = find(mpc.bus(:,2) == PQ);
0120 pv = find(mpc.bus(:,2) == PV);
0121 ref = find(mpc.bus(:,2) == REF);
0122
0123 refpv = sort([ref; pv]);
0124 pvpq = sort([pv; pq]);
0125
0126 Y = makeYbus(mpc);
0127 G = real(Y);
0128 B = imag(Y);
0129
0130
0131 S = makeSbus(mpc.baseMVA, mpc.bus, mpc.gen);
0132 P = real(S);
0133 Q = imag(S);
0134
0135 Vmag = zeros(nbus,1);
0136 Vmag(mpc.gen(:,GEN_BUS)) = mpc.gen(:,VG);
0137
0138
0139
0140
0141
0142 if verbose >= 2
0143 fprintf('Creating variables.\n');
0144 end
0145
0146 Vd = sdpvar(nbus,1);
0147 Vq = sdpvar(nbus,1);
0148
0149 V = [1; Vd; Vq];
0150
0151
0152
0153 if verbose >= 2
0154 fprintf('Creating polynomials.\n');
0155 end
0156
0157 deg = 0;
0158
0159
0160 for i=1:nbus-1
0161 [l,lc] = polynomial(V,deg,0);
0162 lambda(i) = l;
0163 clambda{i} = lc;
0164 end
0165
0166
0167 for i=1:length(pq)
0168 [g,gc] = polynomial(V,deg,0);
0169 gamma(i) = g;
0170 cgamma{i} = gc;
0171 end
0172
0173
0174 for i=1:length(refpv)
0175 [m,mc] = polynomial(V,deg,0);
0176 mu(i) = m;
0177 cmu{i} = mc;
0178 end
0179
0180
0181 [delta, cdelta] = polynomial(V,deg,0);
0182
0183
0184
0185
0186 if verbose >= 2
0187 fprintf('Creating power flow equations.\n');
0188 end
0189
0190
0191 for k=1:nbus-1
0192 i = pvpq(k);
0193
0194
0195
0196
0197 y = find(G(:,i) | B(:,i));
0198 for m=1:length(y)
0199 if exist('fp','var') && length(fp) == k
0200 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)));
0201 else
0202 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)));
0203 end
0204 end
0205
0206 end
0207
0208
0209 for k=1:length(pq)
0210 i = pq(k);
0211
0212
0213
0214
0215 y = find(G(:,i) | B(:,i));
0216 for m=1:length(y)
0217 if exist('fq','var') && length(fq) == k
0218 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)));
0219 else
0220 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)));
0221 end
0222 end
0223 end
0224
0225
0226 for k=1:length(refpv)
0227 i = refpv(k);
0228 fv(k) = Vd(i)^2 + Vq(i)^2;
0229 end
0230
0231
0232
0233
0234 if verbose >= 2
0235 fprintf('Forming the test polynomial.\n');
0236 end
0237
0238 sospoly = 1;
0239
0240
0241 for i=1:nbus-1
0242 sospoly = sospoly + lambda(i) * (fp(i) - P(pvpq(i)));
0243 end
0244
0245
0246 for i=1:length(pq)
0247 sospoly = sospoly + gamma(i) * (fq(i) - Q(pq(i)));
0248 end
0249
0250
0251 for i=1:length(refpv)
0252 sospoly = sospoly + mu(i) * (fv(i) - Vmag(refpv(i))^2);
0253 end
0254
0255
0256 sospoly = sospoly + delta * Vq(ref);
0257
0258 sospoly = -sospoly;
0259
0260
0261
0262 if verbose >= 2
0263 fprintf('Solving the sum of squares problem.\n');
0264 end
0265
0266 sosopts = sdpsettings(sdpopts,'sos.newton',1,'sos.congruence',1);
0267
0268 constraints = sos(sospoly);
0269
0270 pvec = cdelta(:).';
0271
0272 for i=1:length(lambda)
0273 pvec = [pvec clambda{i}(:).'];
0274 end
0275
0276 for i=1:length(gamma)
0277 pvec = [pvec cgamma{i}(:).'];
0278 end
0279
0280 for i=1:length(mu)
0281 pvec = [pvec cmu{i}(:).'];
0282 end
0283
0284
0285 S = warning;
0286 if have_fcn('matlab', 'vnum') >= 8.006 && have_fcn('cplex') && ...
0287 have_fcn('cplex', 'vnum') <= 12.006003
0288 warning('OFF', 'MATLAB:lang:badlyScopedReturnValue');
0289 end
0290
0291 sol = solvesos(constraints,[],sosopts,pvec);
0292
0293 warning(S);
0294
0295 if verbose >= 2
0296 fprintf('Solver exit message: %s\n',sol.info);
0297 end
0298
0299 if sol.problem
0300
0301 insolvable = 0;
0302 elseif sol.problem == 0
0303
0304 insolvable = 1;
0305 end