pinadd2err.java

Go to the documentation of this file.

 
import java.util.*;
import java.lang.Math;
import conopt.*;

public class pinadd2err {
   public static void main(String argv[]) {
      System.loadLibrary("conoptjni4");
 
      String name = "pinadd2err";
 
      Conopt conopt = new Conopt(name);
      MsgHandler msghdlr = new MsgHandler(name);
 
      // adding the message handler to the conopt interface
      conopt.setMessageHandler(msghdlr);
 
      // Turn debugging of FDEval on or off: -1 means first call only
      conopt.debugFV(-1);
 
      // try to set the license using the environment variables
      try {
         int license_int_1 = Integer.parseInt(System.getenv("CONOPT_LICENSE_INT_1"));
         int license_int_2 = Integer.parseInt(System.getenv("CONOPT_LICENSE_INT_2"));
         int license_int_3 = Integer.parseInt(System.getenv("CONOPT_LICENSE_INT_3"));
         String license_text = System.getenv("CONOPT_LICENSE_TEXT");
 
         conopt.setLicense(license_int_1, license_int_2, license_int_3,
               license_text);
      } catch (Exception e) {
         System.out.println("Unable to set license: " + e.getMessage());
      }
 
      std s = new std();
      int retcode = 0;
      int coi_error = 0;
      double[] expectedSols = {1225.8643407378518, 1280.284354018754, 1333.7424507050348, 1386.2354828262864};
 
      for (int i = 16; i <= 20; i++) {
         PinAdd2ErrModelData model = new PinAdd2ErrModelData();
 
         double[] xkeep = conopt.getVariableValues();
         int[] xstat = conopt.getVariableStatus();
         int[] estat = conopt.getConstraintStatus();
 
         // building the model
         model.buildModel(i, xkeep, xstat, estat);
   
         // loading the model in the conopt object
         conopt.loadModel(model);
 
         if (i == 16) {
            conopt.debug2D(-1);
            coi_error = conopt.solve();
 
            if (coi_error != 0) {
               s.java_log(name, "Errors encountered during first solve");
               retcode = -1;
            }
 
            else if (conopt.modelStatus() < 6 || conopt.modelStatus() > 7 ||
                     conopt.solutionStatus() != 5) {
               s.java_log(name, "Incorrect Model or Solver Status during first solve");
               retcode = -1;
            };
 
            conopt.printStatus();
         }
         else {
            /*  Turn debugging of 2nd derivative off again.                       */
            /*  We should expect a proper optimal solution even with an incorrect */
            /*  2DLagr routine so we test for (local) optimality after the solve. */
            conopt.debug2D(0);
 
         conopt.solve();
 
         retcode = s.checkSolve(name, conopt.modelStatus(), conopt.solutionStatus(),
                  conopt.objectiveValue(), expectedSols[i - 17], 0.001);
 
         conopt.printStatus();
         }
      }
      msghdlr.close();
      System.exit(retcode);
   }
}
 
class PinAdd2ErrModelData extends ModelData {
   static final int CONS_SEQ  = 0;
   static final int CONS_DREV = 1;
 
   private int T   = 0;
   private int Vpp = 7;  // Variables per period
   private int Epp = 6;  // Equations per period
 
   // the demand vector
   private double[] demand;
 
   // index arrays for the variables
   private int[] vartd;
   private int[] varcs;
   private int[] vars;
   private int[] vard;
   private int[] varr;
   private int[] varp;
   private int[] varrev;
 
   // index arrays for the constraints
   private int[] consttdeq;
   private int[] constseq;
   private int[] constcseq;
   private int[] constdeq;
   private int[] constreq;
   private int[] constdrev;
 
   private int[] hessianstart;
 
   public class IntPair {
      int first_;
      int second_;
 
      public IntPair(int first, int second) {
         first_ = first;
         second_ = second;
      }
 
      public int first() { return first_; }
      public int second() { return second_; }
   };
 
   private final Map<Integer, IntPair> consmapping = new HashMap<>();
 
 
   public PinAdd2ErrModelData() {
      super();
   }

   public void buildModel(int T, double[] xkeep, int[] xstat, int[] estat) {
 
      this.T = T;
 
      demand = new double[T];
 
      vartd  = new int[T];
      varcs  = new int[T];
      vars   = new int[T];
      vard   = new int[T];
      varr   = new int[T];
      varp   = new int[T];
      varrev = new int[T];
 
      consttdeq = new int[T];
      constseq  = new int[T];
      constcseq = new int[T];
      constdeq  = new int[T];
      constreq  = new int[T];
      constdrev = new int[T];
      hessianstart = new int[T];
 
 
      // adding the variables to the model
      // defining the demand for each time period.
      for (int t = 0; t < T; t++) {
         demand[t] = 1.0 + 2.3 * Math.pow(1.015, t);
      }
 
      // other variables
      int varidx;
      for (int t = 0; t < T; t++) {
         /* td:  Lower=0, Upper=inf, Curr=18 */
         varidx = addVariable(0, Conopt.Infinity, 18);
         vartd[t] = varidx;
 
         /* cs:  Lower=0, Upper=inf, Curr=7*t */
         varidx = addVariable(0, Conopt.Infinity, 7 * (t + 1));
         varcs[t] = varidx;
 
         /* s:   Lower=0, Upper=inf, Curr=7 */
         varidx = addVariable(0, Conopt.Infinity, 7);
         vars[t] = varidx;
 
         /* d:   Lower=0, Upper=inf, Curr=td-s */
         varidx =
               addVariable(0, Conopt.Infinity, getVariable(vartd[t]).getCurr() - getVariable(vars[t]).getCurr());
         vard[t] = varidx;
 
         /* r:   Lower=1, Upper=inf, Curr=r(t-1)-d */
         if (t > 0)
            varidx = addVariable(
                  1, Conopt.Infinity, getVariable(varr[t - 1]).getCurr() - getVariable(vard[t]).getCurr());
         else
            varidx = addVariable(1, Conopt.Infinity, 500 - getVariable(vard[t]).getCurr());
         varr[t] = varidx;
 
         /* p:   Lower=1, Upper=inf, Curr=14 */
         varidx = addVariable(1, Conopt.Infinity, 14);
         varp[t] = varidx;
 
         /* rev: Lower=-inf, Upper=inf, Curr=0 */
         varidx = addVariable(-Conopt.Infinity, Conopt.Infinity);
         varrev[t] = varidx;
      }
 
      if (xkeep.length > 0) {
         for (int i = 0; i < xkeep.length; i++) {
            getVariable(i).setCurr(xkeep[i]);
         }
 
         // linear extrapolation for the last period of variable td
         getVariable(vartd[T - 1]).setCurr(
               2 * getVariable(vartd[T - 2]).getCurr() - getVariable(vartd[T - 3]).getCurr());
         // cs
         getVariable(varcs[T - 1]).setCurr(
               2 * getVariable(varcs[T - 2]).getCurr() - getVariable(varcs[T - 3]).getCurr());
         // s
         getVariable(vars[T - 1]).setCurr(
               2 * getVariable(vars[T - 2]).getCurr() - getVariable(vars[T - 3]).getCurr());
         // d
         getVariable(vard[T - 1]).setCurr(
               2 * getVariable(vard[T - 2]).getCurr() - getVariable(vard[T - 3]).getCurr());
         // r
         getVariable(varr[T - 1]).setCurr(
               2 * getVariable(varr[T - 2]).getCurr() - getVariable(varr[T - 3]).getCurr());
         // p
         getVariable(varp[T - 1]).setCurr(
               2 * getVariable(varp[T - 2]).getCurr() - getVariable(varp[T - 3]).getCurr());
         // rev
         getVariable(varrev[T - 1]).setCurr(
               2 * getVariable(varrev[T - 2]).getCurr() - getVariable(varrev[T - 3]).getCurr());

         for (int i = 0; i < xstat.length; i++) {
            getVariable(i).setVarstatus(xstat[i]);
         }
 
         /* td */
         getVariable(vartd[T - 1]).setVarstatus(getVariable(vartd[T - 2]).getVarstatus());
         /* cs */
         getVariable(varcs[T - 1]).setVarstatus(getVariable(varcs[T - 2]).getVarstatus());
         /* s */
         getVariable(vars[T - 1]).setVarstatus(getVariable(vars[T - 2]).getVarstatus());
         /* d */
         getVariable(vard[T - 1]).setVarstatus(getVariable(vard[T - 2]).getVarstatus());
         /* r */
         getVariable(varr[T - 1]).setVarstatus(getVariable(varr[T - 2]).getVarstatus());
         /* p */
         getVariable(varp[T - 1]).setVarstatus(getVariable(varp[T - 2]).getVarstatus());
         /* rev */
         getVariable(varrev[T - 1]).setVarstatus(getVariable(varrev[T - 2]).getVarstatus());
      }
 
 
      // adding the objective
      double[] objcoeff  = new double[T];
      int[]    objnlflag = new int[T];
 
      for (int t = 0; t < T; ++t) {
         objcoeff[t] = Math.pow(1.05, 1 - (t + 1));
      }
      int objidx = addConstraint(ConstraintType.Free, 0.0, varrev, objcoeff, objnlflag);
 
 
      // adding the constraints to the model
      int considx;
      for (int t = 0; t < T; t++) {
 
         // adding the tdeq equations
         if (t == 0) {
            {
               int[] index = {vartd[t], varp[t]};  // the variables
               double[] value = {1.0, 0.13};       // the coefficients
               int[] nlflag = {0, 0};              // nlflags
               considx = addConstraint(ConstraintType.Eq, demand[t] + 0.87 * 18.0, index, value, nlflag);
            }
         }
         else {
            {
               int[] index = {vartd[t], vartd[t - 1], varp[t]};
               double[] value = {1.0, -0.87, 0.13};
               int[] nlflag = {0, 0, 0};
               considx = addConstraint(ConstraintType.Eq, demand[t], index, value, nlflag);
            }
         }
         consttdeq[t] = considx;
 
         // adding the seq Equations
         if (t == 0) {
            {
               int[] index = {vars[t], varp[t], varcs[t]};
               double[] value = {1.0, 0.0, 0.0};
               int[] nlflag = {0, 1, 1};
               considx = addConstraint(ConstraintType.Eq, 0.75 * 6.5, index, value, nlflag);
            }
         }
         else {
            {
               int[] index = {vars[t], vars[t - 1], varp[t], varcs[t]};
               double[] value = {1.0, -0.75, 0.0, 0.0};
               int[] nlflag = {0, 0, 1, 1};
               considx = addConstraint(ConstraintType.Eq, 0.0, index, value, nlflag);
            }
         }
         constseq[t] = considx;
         consmapping.putIfAbsent(considx, new IntPair(CONS_SEQ, t));
 
         // adding the cseq Equations
         if (t == 0) {
            {
               int[] index = {varcs[t], vars[t]};
               double[] value = {1.0, -1.0};
               int[] nlflag = {0, 0};
               considx = addConstraint(ConstraintType.Eq, 0.0, index, value, nlflag);
            }
         }
         else {
            {
               int[] index = {varcs[t], varcs[t - 1], vars[t]};
               double[] value = {1.0, -1.0, -1.0};
               int[] nlflag = {0, 0, 0};
               considx = addConstraint(ConstraintType.Eq, 0.0, index, value, nlflag);
            }
         }
         constcseq[t] = considx;
 
         // adding the deq equations
         {
               int[] index = {vard[t], vartd[t], vars[t]};
               double[] value = {1.0, -1.0, 1.0};
               int[] nlflag = {0, 0, 0};
               considx = addConstraint(ConstraintType.Eq, 0.0, index, value, nlflag);
         }
         constdeq[t] = considx;
 
         // adding the req equations
         if (t == 0) {
            {
               int[] index = {varr[t], vard[t]};
               double[] value = {1.0, 1.0};
               int[] nlflag = {0, 0};
               considx = addConstraint(ConstraintType.Eq, 500.0, index, value, nlflag);
            }
         }
         else {
            {
               int[] index = {varr[t], varr[t - 1], vard[t]};
               double[] value = {1.0, -1.0, 1.0};
               int[] nlflag = {0, 0, 0};
               considx = addConstraint(ConstraintType.Eq, 0.0, index, value, nlflag);
            }
         }
         constreq[t] = considx;
 
         // adding the drev equations
         {
               int[] index = {varrev[t], vard[t], varp[t], varr[t]};
               double[] value = {1.0, 0.0, 0.0, 0.0};
               int[] nlflag = {0, 1, 1, 1};
               considx = addConstraint(ConstraintType.Eq, 0.0, index, value, nlflag);
         }
         constdrev[t] = considx;
         consmapping.putIfAbsent(considx, new IntPair(CONS_DREV, t));
      }
 
      // Restore values from previous runs
      if (estat.length > 0) {
         // the equation status
         for (int i = 0; i < estat.length; i++) {
            getConstraint(i).setSlackstatus(estat[i]);
         }
 
         // tdeq
         getConstraint(consttdeq[T - 1]).setSlackstatus(getConstraint(consttdeq[T - 2]).getSlackstatus());
         // seq
         getConstraint(constseq[T - 1]).setSlackstatus(getConstraint(constseq[T - 2]).getSlackstatus());
         // cseq
         getConstraint(constcseq[T - 1]).setSlackstatus(getConstraint(constcseq[T - 2]).getSlackstatus());
         // deq
         getConstraint(constdeq[T - 1]).setSlackstatus(getConstraint(constdeq[T - 2]).getSlackstatus());
         // req
         getConstraint(constreq[T - 1]).setSlackstatus(getConstraint(constreq[T - 2]).getSlackstatus());
         // drev
         getConstraint(constdrev[T - 1]).setSlackstatus(getConstraint(constdrev[T - 2]).getSlackstatus());
      }
 
      // setting the objective constraint
      setObjectiveElement(ObjectiveElement.Constraint, objidx);
 
      // setting the optimisation direction
      setOptimizationSense(Sense.Maximize);
 
      // setting the structure of the hessian
      int[] hsrow = new int[T * 5];
      int[] hscol = new int[T * 5];
 
      int idx = 0; // running index into hsrow/hscol
      for (int t = 0; t < T; ++t) {
         hessianstart[t] = idx;
 
         // 1) (cs, cs)
         hscol[idx] = varcs[t];
         hsrow[idx] = varcs[t];
         idx++;
 
         // 2) (cs, p)
         hscol[idx] = varcs[t];
         hsrow[idx] = varp[t];
         idx++;
 
         // 3) (d, r) the row index for this entry is an error.
         hscol[idx] = vard[t];
         hsrow[idx] = varr[t];
         idx++;
 
         // 4) (d, p)
         hscol[idx] = vard[t];
         hsrow[idx] = varp[t];
         idx++;
 
         // 5) (r, r)
         hscol[idx] = varr[t];
         hsrow[idx] = varr[t];
         idx++;
      }
      setSDLagrangianStructure(hsrow, hscol);
   }

   public double evaluateNonlinearTerm(double[] x, int rowno, boolean ignerr, int thread) {
      IntPair conspair = consmapping.get(rowno);
      if (conspair == null) return 0;
 
      int consset = conspair.first();
      int t       = conspair.second();
 
      assert (consset == CONS_SEQ || consset == CONS_DREV);
      double g = 0;
      if (consset == CONS_SEQ) {
         // seq equation. Nonlinear term = -(1.1+0.1*p(t))*1.02**(-cs(t)/7)
         double h1, h2;
 
         h1 = (1.1 + 0.1 * x[varp[t]]);
         h2 = Math.pow(1.02, -x[varcs[t]] / 7.0);
         g = -h1 * h2;
      }
      else if (consset == CONS_DREV) {
         g = -x[vard[t]] * (x[varp[t]] - 250. / x[varr[t]]);
      }
 
      return g;
   }

   public void evaluateNonlinearJacobian(double[] x, double[] jac, int rowno, int[] jacnum, boolean ignerr, int thread) {
      assert x.length == jac.length;
 
      IntPair conspair = consmapping.get(rowno);
      if (conspair == null) return;
 
      int consset = conspair.first();
      int t       = conspair.second();
 
      assert (consset == CONS_SEQ || consset == CONS_DREV);
 
      if (consset == CONS_SEQ) {
         // seq equation. Nonlinear term = -(1.1+0.1*p(t))*1.02**(-cs(t)/7)
         double h1, h2;
 
         h1 = (1.1 + 0.1 * x[varp[t]]);
         h2 = Math.pow(1.02, -x[varcs[t]] / 7.0);
 
         jac[varcs[t]] = h1 * h2 * Math.log(1.02) / 7.0;
         jac[varp[t]] = -h2 * 0.1;
      }
      else if (consset == CONS_DREV) {
         jac[vard[t]] = -(x[varp[t]] - 250. / x[varr[t]]);
         jac[varr[t]] = -x[vard[t]] * 250. / Math.pow(x[varr[t]], 2);
         jac[varp[t]] = -x[vard[t]];
      }
   }

   public void evaluateSDLagrangian(double x[], double u[], int[] hessianrow, int[] hessiancol, double[] hessianval) {
 
      double cs, d, r, p;
      double h1, h2;
 
      /*  Normal Evaluation mode                                            */
      for (int t = 0; t < T; t++) {
         int hessidx = hessianstart[t];
         cs = x[varcs[t]];
         d = x[vard[t]];
         r = x[varr[t]];
         p = x[varp[t]];
 
         h1 = 1.1 + 0.1 * p;
         h2 = Math.pow(1.02, -cs / 7.0);
 
         /* the second derivative of the SEQ equations */
         hessianval[hessidx] = -h1 * h2 * Math.pow(Math.log(1.02) / 7.0, 2) * u[constseq[t]];
         hessianval[hessidx + 1] = 0.1 * h2 * (Math.log(1.02) / 7.) * u[constseq[t]];
 
         /* the second derivative of the DREV equations */
         hessianval[hessidx + 2] = -250.0 / Math.pow(r, 2) * u[constdrev[t]];
         hessianval[hessidx + 3] = -1.0 * u[constdrev[t]];
         hessianval[hessidx + 4] = 500.0 * d / Math.pow(r, 3) * u[constdrev[t]];
      }
   }
}