Code conversion from GAMS to JuMP (Ipopt Solver)

SingingKim · July 31, 2017, 1:43pm

Is there any better way to convert constraints with multiple indices?
Here’s what I’ve done.

[GAMS]
SETS
         i       index of buses         /1*33/
         t       index of time slots    /1*24/
         s      reduced scenarios      /1*10/
;

* Power flow equations
* Bus: 2~N_Bus
NLECONS1(t,i,s)$(ord(i) > 1) .. p_g(t,i,s)+p_bes(t,i,s)-p_l(t,i,s)-v(t,i,s)*SUM(j, v(t,j,s)*(GmatDATA(i,j)*cos((2*pi/360)*(theta(t,i,s)-theta(t,j,s)))+BmatDATA(i,j)*sin((2*pi/360)*(theta(t,i,s)-theta(t,j,s)))))=e=0   ;
NLECONS2(t,i,s)$(ord(i) > 1) .. q_g(t,i,s)+q_bes(t,i,s)-q_l(t,i,s)-v(t,i,s)*SUM(j, v(t,j,s)*(GmatDATA(i,j)*sin((2*pi/360)*(theta(t,i,s)-theta(t,j,s)))-BmatDATA(i,j)*cos((2*pi/360)*(theta(t,i,s)-theta(t,j,s)))))=e=0   ;
* Slack Bus
NLECONS3(t,i,s)$(ord(i) = 1) .. p_tie(t,s)+p_g(t,i,s)+p_bes(t,i,s)-p_l(t,i,s)-v(t,i,s)*SUM(j, v(t,j,s)*(GmatDATA(i,j)*cos((2*pi/360)*(theta(t,i,s)-theta(t,j,s)))+BmatDATA(i,j)*sin((2*pi/360)*(theta(t,i,s)-theta(t,j,s)))))=e=0  ;
NLECONS4(t,i,s)$(ord(i) = 1) .. q_tie(t,s)+q_g(t,i,s)+q_bes(t,i,s)-q_l(t,i,s)-v(t,i,s)*SUM(j, v(t,j,s)*(GmatDATA(i,j)*sin((2*pi/360)*(theta(t,i,s)-theta(t,j,s)))-BmatDATA(i,j)*cos((2*pi/360)*(theta(t,i,s)-theta(t,j,s)))))=e=0  ;

[JuMP]
# Bus: 2~N_Bus
ConstraintIndex4 = Array{Tuple{Int64,Int64,Int64}}(N_Scen,N_TimeSlot,N_Bus-1);
for i in 1:N_Scen
  for j in 1:N_TimeSlot
    for k in 2:N_Bus
      ConstraintIndex4[i,j,k-1] = (i,j,k)
    end
  end
end
@NLconstraint( m, POWERFLOWconstraint1[ (s,t,i) in ConstraintIndex4], Pgen[s,t,i] +Pbess[s,t,i] -Pload[s,t,i] -v[s,t,i] * sum(v[s,t,j]*(data["Gmat"][i,j]*cos((2*pi/360)*(theta[s,t,i]-theta[s,t,j]))+data["Bmat"][i,j]*sin((2*pi/360)*(theta[s,t,i]-theta[s,t,j]))) for j in 1:N_Bus)==0)
@NLconstraint( m, POWERFLOWconstraint2[ (s,t,i) in ConstraintIndex4], Qgen[s,t,i] +Qbess[s,t,i] -Qload[s,t,i] -v[s,t,i] * sum(v[s,t,j]*(data["Gmat"][i,j]*sin((2*pi/360)*(theta[s,t,i]-theta[s,t,j]))-data["Bmat"][i,j]*cos((2*pi/360)*(theta[s,t,i]-theta[s,t,j]))) for j in 1:N_Bus)==0)

# Slack BUS
ConstraintIndex = Array{Tuple{Int64,Int64}}(N_Scen,N_TimeSlot);
for i in 1:N_Scen
  for j in 1:N_TimeSlot
    ConstraintIndex[i,j] = (i,j)
  end
end
@NLconstraint( m, POWERFLOWconstraint3[ (s,t) in ConstraintIndex], Ptie[s,t]+Pgen[s,t,1] +Pbess[s,t,1] -Pload[s,t,1] -v[s,t,1] * sum(v[s,t,j]*(data["Gmat"][1,j]*cos((2*pi/360)*(theta[s,t,1]-theta[s,t,j]))+data["Bmat"][1,j]*sin((2*pi/360)*(theta[s,t,1]-theta[s,t,j]))) for j in 1:N_Bus)==0)
@NLconstraint( m, POWERFLOWconstraint4[ (s,t) in ConstraintIndex], Qtie[s,t]+Qgen[s,t,1] +Qbess[s,t,1] -Qload[s,t,1] -v[s,t,1] * sum(v[s,t,j]*(data["Gmat"][1,j]*sin((2*pi/360)*(theta[s,t,1]-theta[s,t,j]))-data["Bmat"][1,j]*cos((2*pi/360)*(theta[s,t,1]-theta[s,t,j]))) for j in 1:N_Bus)==0)

Ipopt Solvers on GAMS and JuMP show different numbers of nonzeros in equality constraint Jacobian even though total numbers of equality constraints are same.


[GAMS]
Number of nonzeros in equality constraint Jacobian...:   303010  
Total number of variables............................:    66864
Total number of equality constraints.................:    39990
Total number of inequality constraints...............:    32640

[JuMP]
Number of nonzeros in equality constraint Jacobian...:  1245790
Total number of variables............................:    66864
Total number of equality constraints.................:    39990
Total number of inequality constraints...............:    32640

GAMS solves the problem in 1~2hrs, but JuMP couldn’t finish the work in days.
I don’t know where this comes from.

[GAMS]

COIN-OR Interior Point Optimizer (Ipopt Library 3.8)
written by A. Waechter

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Common Public License (CPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

NOTE: You are using Ipopt by default with the MUMPS linear solver.
      Other linear solvers might be more efficient (see Ipopt documentation).


This is Ipopt version 3.8stable, running with linear solver mumps.

Number of nonzeros in equality constraint Jacobian...:   303010
Number of nonzeros in inequality constraint Jacobian.:    47100
Number of nonzeros in Lagrangian Hessian.............:    87794

Total number of variables............................:    66864
                     variables with only lower bounds:     8480
                variables with lower and upper bounds:    34100
                     variables with only upper bounds:        0
Total number of equality constraints.................:    39990
Total number of inequality constraints...............:    32640
        inequality constraints with only lower bounds:     8640
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:    24000

[JuMP]

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.4, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:  1245790
Number of nonzeros in inequality constraint Jacobian.:    72560
Number of nonzeros in Lagrangian Hessian.............: 35581600

Total number of variables............................:    66864
                     variables with only lower bounds:     8480
                variables with lower and upper bounds:    34100
                     variables with only upper bounds:        0
Total number of equality constraints.................:    39990
Total number of inequality constraints...............:    32640
        inequality constraints with only lower bounds:     8640
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:    24000

tkelman · August 1, 2017, 9:39am

JuMP doesn’t do much preprocessing to combine duplicate elements, it makes the solver do some of the work there. Looks like the Hessian in particular is handled better in GAMS than JuMP. You could try the AmplNLWriter package.

SingingKim · August 1, 2017, 3:49pm

Thanks for the comment.
Doesn’t AmplNLWriter calls Ipopt Solver for NLP?
As far as I recall, BonminSolver does.

tkelman · August 2, 2017, 5:18am

It’s the same optimization solver but uses a C library for evaluating the Jacobian, Hessian, etc instead of the Julia code used by Ipopt.jl

SingingKim · August 2, 2017, 7:07am

@tkelman Thanks for the kind comments and advice.
I ran the same code with AmplNLSolver("Ipopt") .
It’s been 2hrs and still running (GAMS Ipopt completes the computation in 1~2hrs), and giving no variable/constraint summary but some expressions deprecation warnings.
Doesn’t AmplNLWriter package provide iteration status?

tkelman · August 2, 2017, 12:49pm

Each solver in Julia handles options slightly differently. There is a way to turn on solver logging at intermediate iterations but it may be spelled a little differently in each.