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