src/main/java/edu/cornell/pserc/jpower/opf/Djp_jipsopf_solver.java

   1 /*
   2  * Copyright (C) 1996-2010 Power System Engineering Research Center
   3  * Copyright (C) 2010-2011 Richard Lincoln
   4  *
   5  * JPOWER is free software: you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published
   7  * by the Free Software Foundation, either version 3 of the License,
   8  * or (at your option) any later version.
   9  *
  10  * JPOWER is distributed in the hope that it will be useful,
  11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  13  * GNU General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with JPOWER. If not, see <http://www.gnu.org/licenses/>.
  17  *
  18  */
  19
  20 package edu.cornell.pserc.jpower.opf;
  21
  22 import java.util.HashMap;
  23 import java.util.Map;
  24
  25 import cern.colt.matrix.AbstractMatrix;
  26 import cern.colt.matrix.tdcomplex.DComplexMatrix1D;
  27 import cern.colt.matrix.tdcomplex.DComplexMatrix2D;
  28 import cern.colt.matrix.tdouble.DoubleFactory1D;
  29 import cern.colt.matrix.tdouble.DoubleFactory2D;
  30 import cern.colt.matrix.tdouble.DoubleMatrix1D;
  31 import cern.colt.matrix.tdouble.DoubleMatrix2D;
  32
  33 import static edu.emory.mathcs.utils.Utils.ifunc;
  34 import static edu.emory.mathcs.utils.Utils.dfunc;
  35 import static edu.emory.mathcs.utils.Utils.cfunc;
  36 import static edu.emory.mathcs.utils.Utils.intm;
  37 import static edu.emory.mathcs.utils.Utils.irange;
  38 import static edu.emory.mathcs.utils.Utils.nonzero;
  39 import static edu.emory.mathcs.utils.Utils.complex;
  40
  41 import static edu.cornell.pserc.jpower.jpc.JPC.PW_LINEAR;
  42 import static edu.cornell.pserc.jpower.jpc.JPC.REF;
  43 import static edu.cornell.pserc.jpower.pf.Djp_makeYbus.makeYbus;
  44 import static edu.cornell.pserc.jips.Dips_jips.jips;
  45
  46 import edu.cornell.pserc.jips.ConstraintEvaluator;
  47 import edu.cornell.pserc.jips.HessianEvaluator;
  48 import edu.cornell.pserc.jips.ObjectiveEvaluator;
  49 import edu.cornell.pserc.jpower.jpc.Branch;
  50 import edu.cornell.pserc.jpower.jpc.Bus;
  51 import edu.cornell.pserc.jpower.jpc.Gen;
  52 import edu.cornell.pserc.jpower.jpc.Cost;
  53 import edu.cornell.pserc.jpower.jpc.JPC;
  54 import edu.cornell.pserc.jpower.opf.Djp_opf_model.Set;
  55
  56 /**
  57  * Solves AC optimal power flow using JIPS.
  58  *
  59  * @author Ray Zimmerman
  60  * @author Richard Lincoln
  61  *
  62  */
  63 public class Djp_jipsopf_solver {
  64
  65         /**
  66          * Solves AC optimal power flow using JIPS.
  67          *
  68          * @param om
  69          * @param jpopt
  70          * @param output
  71          * @return
  72          */
  73         @SuppressWarnings("static-access")
  74         public static Object[] jipsopf_solver(Djp_opf_model om,
  75                         Map<String, Double> jpopt, Map<String, AbstractMatrix> out_opt) {
  76
  77                 int nb, nl, ny, nl2, info, nlnN;
  78                 int[] ipwl, il, kl, ku;
  79                 boolean verbose, success;
  80                 double feastol, gradtol, comptol, costtol, max_it, max_red,
  81                                 step_control, baseMVA, c, f;
  82                 Map<String, Double> opt;
  83                 JPC jpc, results;
  84                 Bus bus;
  85                 Gen gen;
  86                 Branch branch;
  87                 Cost gencost;
  88
  89                 Map<String, Set>[] idx;
  90                 Map<String, Set> vv, ll, nn;
  91                 Map<String, Object> Output, raw;
  92                 Map<String, DoubleMatrix1D> Lambda, var, nln, lin;
  93                 Map<String, Map<String, DoubleMatrix1D>> mu;
  94
  95                 AbstractMatrix[] Alu;
  96                 DoubleMatrix1D l, u, x0, xmin, xmax, lb, ub, Varefs, x,
  97                                 Va, Vm, Pg, Qg, muSf, muSt, nl_mu_l, nl_mu_u, pimul;
  98                 DoubleMatrix1D[] xx;
  99                 DoubleMatrix1D[] geq;
 100                 DoubleMatrix2D A;
 101                 DComplexMatrix1D V, Sf, St;
 102                 DComplexMatrix2D Ybus, Yf, Yt;
 103                 DComplexMatrix2D[] Y;
 104
 105                 ObjectiveEvaluator f_fcn;
 106                 HessianEvaluator hess_fcn;
 107                 ConstraintEvaluator gh_fcn, geq_fcn;
 108
 109                 Object[] jips;
 110
 111                 /* ----- initialization ----- */
 112
 113                 /* options */
 114                 verbose = jpopt.get("VERBOSE") == 1;
 115                 feastol = jpopt.get("PDIPM_FEASTOL");
 116                 gradtol = jpopt.get("PDIPM_GRADTOL");
 117                 comptol = jpopt.get("PDIPM_COMPTOL");
 118                 costtol = jpopt.get("PDIPM_COSTTOL");
 119                 max_it  = jpopt.get("PDIPM_MAX_IT");
 120                 max_red = jpopt.get("SCPDIPM_RED_IT");
 121                 step_control = jpopt.get("OPF_ALG") == 565 ? 1:0;
 122                 if (feastol == 0)
 123                         feastol = jpopt.get("OPF_VIOLATION");   // = OPF_VIOLATION by default
 124                 opt = new HashMap<String, Double>();
 125                 opt.put("feastol", feastol);
 126                 opt.put("gradtol", gradtol);
 127                 opt.put("comptol", comptol);
 128                 opt.put("costtol", costtol);
 129                 opt.put("max_it", max_it);
 130                 opt.put("max_red", max_red);
 131                 opt.put("step_control", step_control);
 132                 opt.put("cost_mult", 1e-4);
 133                 opt.put("verbose", verbose ? 1.0 : 0.0);
 134
 135                 /* unpack data */
 136                 jpc = om.get_jpc();
 137                 baseMVA = jpc.baseMVA;
 138                 bus = jpc.bus;
 139                 gen = jpc.gen;
 140                 branch = jpc.branch;
 141                 gencost = jpc.gencost;
 142                 idx = om.get_idx();
 143                 vv = idx[0]; ll = idx[1]; nn = idx[2];
 144
 145                 /* problem dimensions */
 146                 nb = bus.size();                        // number of buses
 147                 nl = branch.size();                     // number of branches
 148                 ny = om.getN("var", "y");       // number of piece-wise linear costs
 149
 150                 /* linear constraints */
 151                 Alu = om.linear_constraints();
 152                 A = (DoubleMatrix2D) Alu[0];
 153                 l = (DoubleMatrix1D) Alu[1];
 154                 u = (DoubleMatrix1D) Alu[2];
 155
 156                 xx = om.getv();
 157                 x0 = xx[0]; xmin = xx[1]; xmax = xx[2];
 158
 159                 /* build admittance matrices */
 160                 Y = makeYbus(baseMVA, bus, branch);
 161                 Ybus = Y[0]; Yf = Y[1]; Yt = Y[2];
 162
 163                 /* try to select an interior initial point */
 164                 lb = xmin.copy(); ub = xmax.copy();
 165                 lb.viewSelection( intm( xmin.copy().assign(dfunc.equals(Double.NEGATIVE_INFINITY)) ).toArray() ).assign(-1e10); // replace Inf with numerical proxies
 166                 ub.viewSelection( intm( xmax.copy().assign(dfunc.equals(Double.POSITIVE_INFINITY)) ).toArray() ).assign( 1e10);
 167                 x0 = lb.copy().assign(ub, dfunc.plus).assign(dfunc.div(2));
 168                 Varefs = bus.Va.viewSelection( bus.bus_type.copy().assign(ifunc.equals(REF)).toArray() ).assign(dfunc.mult(Math.PI)).assign(dfunc.div(180));
 169                 x0.viewSelection(irange(vv.get("Va").i0, vv.get("Va").iN)).assign(Varefs.get(0));       // angles set to first reference angle
 170                 if (ny > 0) {
 171                         ipwl = nonzero( gencost.model.copy().assign(ifunc.equals(PW_LINEAR)) );
 172                         c = gencost.cost.viewSelection(ipwl, null).getMaxLocation()[0];         // largest y-value in CCV data
 173                         // TODO: compute c using sub2ind
 174                         x0.viewPart(vv.get("y").i0, vv.get("y").N).assign(c + 0.1 * Math.abs(c));
 175                 }
 176
 177                 /* find branches with flow limits */
 178                 il = nonzero( intm( branch.rate_a.copy().assign(dfunc.equals(0)) ).assign(ifunc.not).assign(intm( branch.rate_a.assign(dfunc.less(1e10)) ), ifunc.and) );
 179                 nl2 = il.length;        // number of constrained lines
 180
 181                 /* -----  run opf  ----- */
 182                 f_fcn = new Djp_opf_costfcn(om);
 183                 hess_fcn = new Djp_opf_hessfcn(om, Ybus, Yf.viewSelection(il, null), Yt.viewSelection(il, null), jpopt, il, opt.get("cost_mult"));
 184                 gh_fcn = new Djp_opf_consfcn(om, Ybus, Yf.viewSelection(il, null), Yt.viewSelection(il, null), jpopt, il);
 185                 jips = jips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
 186                 x = (DoubleMatrix1D) jips[0];
 187                 f = (Double) jips[1];
 188                 info = (Integer) jips[2];
 189                 Output = (Map<String, Object>) jips[3];
 190                 Lambda = (Map<String, DoubleMatrix1D>) jips[4];
 191
 192                 success = (info > 0);
 193
 194                 /* update solution data */
 195                 Va = x.viewPart(vv.get("Va").i0, vv.get("Va").N).copy();
 196                 Vm = x.viewPart(vv.get("Vm").i0, vv.get("Vm").N).copy();
 197                 Pg = x.viewPart(vv.get("Pg").i0, vv.get("Pg").N).copy();
 198                 Qg = x.viewPart(vv.get("Qg").i0, vv.get("Qg").N).copy();
 199                 V  = complex(Vm, Va);
 200
 201                 /* -----  calculate return values  ----- */
 202
 203                 /* update voltages & generator outputs */
 204                 bus.Va = Va.assign(dfunc.mult(180)).assign(dfunc.div(Math.PI));
 205                 bus.Vm = Vm;
 206                 gen.Pg = Pg.assign(dfunc.mult(baseMVA));
 207                 gen.Qg = Qg.assign(dfunc.mult(baseMVA));
 208                 gen.Vg = Vm.viewSelection(gen.gen_bus.toArray());
 209
 210                 /* compute branch flows */
 211                 Sf = Yf.zMult(V, null).assign(cfunc.conj).assign(V.viewSelection(branch.f_bus.toArray()), cfunc.mult);  // cplx pwr at "from" bus, p.u.
 212                 St = Yt.zMult(V, null).assign(cfunc.conj).assign(V.viewSelection(branch.t_bus.toArray()), cfunc.mult);  // cplx pwr at "to" bus, p.u.
 213                 branch.Pf = Sf.getRealPart().assign(dfunc.mult(baseMVA));
 214                 branch.Qf = Sf.getImaginaryPart().assign(dfunc.mult(baseMVA));
 215                 branch.Pt = St.getRealPart().assign(dfunc.mult(baseMVA));
 216                 branch.Qt = St.getImaginaryPart().assign(dfunc.mult(baseMVA));
 217
 218                 /* line constraint is actually on square of limit, so we must fix multipliers */
 219                 muSf = DoubleFactory1D.dense.make(nl);
 220                 muSt = DoubleFactory1D.dense.make(nl);
 221                 if (il.length > 0) {
 222                         muSf.viewSelection(il).assign( Lambda.get("ineqnonlin").viewPart(0,   nl2).copy().assign(branch.rate_a.viewSelection(il), dfunc.mult).assign(dfunc.mult(2)).assign(dfunc.div(baseMVA)) );
 223                         muSt.viewSelection(il).assign( Lambda.get("ineqnonlin").viewPart(nl2, nl2).copy().assign(branch.rate_a.viewSelection(il), dfunc.mult).assign(dfunc.mult(2)).assign(dfunc.div(baseMVA)) );
 224                 }
 225
 226                 /* update Lagrange multipliers */
 227                 bus.mu_Vmax = Lambda.get("upper").viewPart(vv.get("Vm").i0, vv.get("Vm").N).copy();
 228                 bus.mu_Vmin = Lambda.get("lower").viewPart(vv.get("Vm").i0, vv.get("Vm").N).copy();
 229                 gen.mu_Pmax = Lambda.get("upper").viewPart(vv.get("Pg").i0, vv.get("Pg").N).copy().assign(dfunc.div(baseMVA));
 230                 gen.mu_Pmin = Lambda.get("lower").viewPart(vv.get("Pg").i0, vv.get("Pg").N).copy().assign(dfunc.div(baseMVA));
 231                 gen.mu_Qmax = Lambda.get("upper").viewPart(vv.get("Qg").i0, vv.get("Qg").N).copy().assign(dfunc.div(baseMVA));
 232                 gen.mu_Qmin = Lambda.get("lower").viewPart(vv.get("Qg").i0, vv.get("Qg").N).copy().assign(dfunc.div(baseMVA));
 233                 bus.lam_P = Lambda.get("eqnonlin").viewPart(nn.get("Pmis").i0, nn.get("Pmis").N).copy().assign(dfunc.div(baseMVA));
 234                 bus.lam_Q = Lambda.get("eqnonlin").viewPart(nn.get("Qmis").i0, nn.get("Qmis").N).copy().assign(dfunc.div(baseMVA));
 235                 branch.mu_Sf = muSf.copy().assign(dfunc.div(baseMVA));
 236                 branch.mu_St = muSt.copy().assign(dfunc.div(baseMVA));
 237
 238                 /* package up results */
 239                 nlnN = om.getN("nln");
 240
 241                 /* extract multipliers for non-linear constraints */
 242                 kl = nonzero( Lambda.get("eqnonlin").copy().assign(dfunc.less(0)) );
 243                 ku = nonzero( Lambda.get("eqnonlin").copy().assign(dfunc.greater(0)) );
 244                 nl_mu_l = DoubleFactory1D.dense.make(nlnN);
 245                 nl_mu_u = DoubleFactory1D.dense.make(new DoubleMatrix1D[] {DoubleFactory1D.dense.make(nb), muSf, muSt});
 246                 nl_mu_l.viewSelection(kl).assign( Lambda.get("eqnonlin").viewSelection(kl).copy().assign(dfunc.neg) );
 247                 nl_mu_u.viewSelection(ku).assign( Lambda.get("eqnonlin").viewSelection(ku) );
 248
 249                 mu = new HashMap<String, Map<String,DoubleMatrix1D>>();
 250                 var = new HashMap<String, DoubleMatrix1D>();
 251                 var.put("l", Lambda.get("lower"));
 252                 var.put("u", Lambda.get("upper"));
 253                 mu.put("var", var);
 254                 nln = new HashMap<String, DoubleMatrix1D>();
 255                 nln.put("l", nl_mu_l);
 256                 nln.put("u", nl_mu_u);
 257                 mu.put("nln", nln);
 258                 lin = new HashMap<String, DoubleMatrix1D>();
 259                 lin.put("l", Lambda.get("mu_l"));
 260                 lin.put("u", Lambda.get("mu_u"));
 261                 mu.put("lin", lin);
 262
 263                 results = jpc.copy();
 264                 results.bus = bus.copy();
 265                 results.branch = branch.copy();
 266                 results.gen = gen.copy();
 267                 results.om = om;
 268                 results.x = x.copy();
 269                 results.mu = mu;
 270                 results.f = f;
 271
 272                 /* optional fields */
 273                 if (out_opt.containsKey("dg")) {
 274                         geq_fcn = new Djp_opf_consfcn(om, Ybus, Yf, Yt, jpopt);
 275                         geq = geq_fcn.gh(x);
 276                         results.g  = DoubleFactory1D.dense.append(geq[1], geq[0]);      // include this since we computed it anyway
 277                         DoubleMatrix2D[] dgeq = geq_fcn.dgh(x);
 278                         results.dg = DoubleFactory2D.sparse.appendRows(dgeq[1].viewDice(), dgeq[0].viewDice()); // true Jacobian organization
 279                 }
 280                 if (out_opt.containsKey("g") && !out_opt.containsKey("dg")) {
 281                         geq_fcn = new Djp_opf_consfcn(om, Ybus, Yf, Yt, jpopt);
 282                         geq = geq_fcn.gh(x);
 283                         results.g = DoubleFactory1D.dense.append(geq[1], geq[0]);
 284                 }
 285                 if (out_opt.containsKey("df"))
 286                         results.df = null;
 287                 if (out_opt.containsKey("d2f"))
 288                         results.d2f = null;
 289
 290                 pimul = DoubleFactory1D.dense.make(new DoubleMatrix1D[] {
 291                                 results.mu.get("nln").get("l").copy().assign(results.mu.get("nln").get("u"), dfunc.minus),
 292                                 results.mu.get("lin").get("l").copy().assign(results.mu.get("lin").get("u"), dfunc.minus),
 293                                 DoubleFactory1D.dense.make((ny>0) ? 1:0, -1),
 294                                 results.mu.get("var").get("l").copy().assign(results.mu.get("var").get("u"), dfunc.minus)
 295                 });
 296
 297                 raw = new HashMap<String, Object>();
 298                 raw.put("xr", x);
 299                 raw.put("pimul", pimul);
 300                 raw.put("info", info);
 301                 raw.put("output", Output);
 302
 303                 return new Object[] {results, success, raw};
 304         }
 305
 306         public static Object[] jipsopf_solver(Djp_opf_model om, Map<String, Double> jpopt) {
 307                 return jipsopf_solver(om, jpopt, new HashMap<String, AbstractMatrix>());
 308         }
 309
 310 }