Problems Solving Quadratic Equation Using Ipopt

[BA]

Hi All,

Im am working through an optimization problem which i mentioned a few days ago.  I got everything to run in ipopt but the solver is giving me some perplexing results (again im a newbie so Im sure im making some obvious mistakes).  

The objective is to maximize revenue as defined as:

Pricew*Demandw + Pricesalvage*Inventoryfinal.

I want to optimize the former term on weekly basis and iteratively such that the optimal price for week,w, factors in the subsequent price/demand dynamics for subsequent weeks.  The summary the mathematical objective is:

Max  ΣPwλw + ΣPsalvgInv_end

Where: 

• Pw = price in week, w;  

• λw = demand in week, w;

• Psalv = Salvage price

• Inv_end = Final salvage inventory

Here is my model thus far:

model = AbstractModel()

*******Params and Vars *******************

# Number of weeks

model.Wks = RangeSet(1,12)  

# Salvage price

model.Psalv = Param(initialize=15)  

# Inventory at end of week, w

model.Inv_end = Var(model.Wks, within=NonNegativeIntegers, bounds=(0,2000), initialize=inv)

# Demand in week, w  

model.Wk_dmd = Var(model.Wks, within=NonNegativeIntegers, bounds=(0,2000)) 

#Price in week , w

model.Pwk = Var(model.Wks, within=NonNegativeReals, initialize=5)

*******Objective*******************

def objective_rule(model):

    return sum((model.Wk_dmd[w]*model.Pwk[w] for w in model.Wks) + model.Psalv * model.Inv_end[12])

model.objective = Objective(rule=objective_rule, sense=maximize)

************Constraints**********************

def Dmd_dynamics_rule(model, w):

   ###Calcs demand for given week###

   if w == 1:

      return model.Wk_dmd[w] == inv*.15          

   else:

       # Main demand equation

      return model.Wk_dmd[w] == model.Wk_dmd[w-1]*.3+og_ct*.001+30*(200/age)+.01*inv+2000/model.Pwk[w]   #vars og_ct, age, and inv are global vats

model.Dmd_dynamics = Constraint(model.Wks,rule=Dmd_dynamics_rule)

def Inv_dmd_dynam_rule(model, w):

    ### Prev w End inv is ceiling to current w demand###

    if w > 1:

        return (model.Wk_dmd[w] <= model.Inv_end[w-1])

    return Constraint.Skip

model.Inv_dmd_dynam =Constraint(model.Wks, rule=Inv_dmd_dynam_rule)                                

def Price_dynam_rule(model, w):

     ###Price can not go up, only down###

     if w > 1:

          return model.Pcrrnt[w] <= model.Pcrrnt[w-1]

     return Constraint.Skip

model.Price_dynam = Constraint(model.Wks,rule=Price_dynam_rule)

def Calc_inv(model, w):

    ### Calculates Ending inventory in week, w ###

    if w > 1:

          return model.Inv_end[w] == model.Inv_end[w-1] - model.Wk_dmd[w]

    return Constraint.Skip

model.Inv_end_dynamics = Constraint(model.Wks, rule=Calc_inv).

A couple takeaways from the solver results that I thought were peculiar:

• Objective is listed as 'None' in the suffix declarations
• 'sense' is showing as 'unknown' though i specify sense=maximize in the objective function
• The solution appears to stop after only 1 iteration through the vars specifying "Maximum Number of Iterations Exceeded".  How do I increase this.
• My Inv_End value starts at 7.99673555664e+13 despite an upper bound of 2,000

Here is some additional detail on the solution output below:

 = Solver Results                                         =

# ==========================================================

# ----------------------------------------------------------

#   Problem Information

# ----------------------------------------------------------

Problem: 

- Lower bound: -inf

  Upper bound: inf

  Number of objectives: 1

  Number of constraints: 21

  Number of variables: 18

  Sense: unknown

# ----------------------------------------------------------

#   Solver Information

# ----------------------------------------------------------

Solver: 

- Status: warning

  Message: Ipopt 3.11.9\x3a Maximum Number of Iterations Exceeded.

  Termination condition: maxIterations

  Id: 400

  Error rc: 0

# ----------------------------------------------------------

#   Solution Information

# ----------------------------------------------------------

Solution: 

- number of solutions: 1

  number of solutions displayed: 1

- Status: stoppedByLimit

  Message: Ipopt 3.11.9\x3a Maximum Number of Iterations Exceeded.

  Variable: 

    v0:

      Id: 0

      Value: 7.99673555664e+13

    v1:

      Id: 1

      Value: 6.90382407242e+13

    v2:

      Id: 2

      Value: 6.08140328328e+13

    v3:

      Id: 3

      Value: 5.42863049624e+13

    v4:

      Id: 4

      Value: 4.92201920122e+13

    v5:

      Id: 5

      Value: 1.75848644786e+14

    v6:

      Id: 6

      Value: 409.812936602

    v7:

      Id: 7

      Value: 264.736414547

    v8:

      Id: 8

      Value: 204.692700811

    v9:

      Id: 9

      Value: 175.147505697

    v10:

      Id: 10

      Value: 156.440929811

    v11:

      Id: 11

      Value: 136.42489406

    v12:

      Id: 12

      Value: 1999.99983144

    v13:

      Id: 13

      Value: 1735.26341689

    v14:

      Id: 14

      Value: 1530.57071608

    v15:

      Id: 15

      Value: 1355.42321038

    v16:

      Id: 16

      Value: 1198.98228057

    v17:

      Id: 17

      Value: 1062.55738651

  Constraint: No nonzero values

  Status description: EXCEEDED MAXIMUM NUMBER OF ITERATIONS\x3a the solver was stopped by a limit that you set!

[Watson, Jean-Paul]

You should try to solve with the option “—stream-solver”, assuming you are running with the “pyomo” command-line tool. Basically, it looks like ipopt is not converging – this is a model/ipopt issue, not a pyomo issue.

jpw

--
You received this message because you are subscribed to the Google Groups "Pyomo Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Blake Adams]

Thanks for the reply.  So my model construction is sound?

Also I'm actually using pyomo scripting - is there an equivalent means to get the stream-solver option?

You received this message because you are subscribed to a topic in the Google Groups "Pyomo Forum" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/pyomo-forum/ggvPaZ1HTk4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to pyomo-forum...@googlegroups.com.

[Carl Laird]

Hi,

You can turn on the solver streaming with the option “tee=True” in the opt.solve call (see below)

opt = SolverFactory('ipopt',solver_io='nl')
results = opt.solve(instance,
                    tee=True)

This will show you the output from ipopt.

Regards,

Carl D. Laird
Associate Professor, School of Chemical Engineering, Purdue University
Faculty Fellow, Mary Kay O’Connor Process Safety Center

Ph: 765.494.0085
Email: carl...@purdue.edu

[Blake Adams]

Thanks you  - very much appreciated.

[Watson, Jean-Paul]

I can’t guarantee that it’s sound (in the sense that the equations are all correct), but it can’t be horribly wrong – ipopt understood and parsed what was written.

Because you are scripting, it may be a good idea to invoke “pprint()” on the model instance after it is constructed, and prior to the solve. Assuming it isn’t too large, you should be able to validate that the index sets, variables, and constraints that it reports – in a somewhat human readable form – are what you intended. 

Beyond that, model debugging and validation can be a rather manual and intensive process. You may still have some odd issue with your model, or it could be ipopt. More specifically, it may be the starting point that ipopt is using. Which, because you aren’t supplying one (I think), defaults to something like 0.0 or the mid-point between variable slower and upper bounds. You don’t have to care about the initial point for linear programs, but you do for nonlinear solvers. 

As for the option you need to supply to the solve() command, it is “tee=True”.

jpw