Choice variable doesn't optimize the objective function

Ayesha · February 23, 2024, 10:21am

begin # set parameters
	b = 100.0
	r = 0.2 
	α = 0.85
	θ = 2.5
	β = 0.9
	A = 1.0 
	s = 1.0 
end

begin 
	f(b) = b^α
	Y₂(b,θ) = θ*f(b)
	C₂(R₂,b,θ) = Y₂(b,θ) - R₂
	R̄(b,r) = (1+r)*b
	U₂(R₂, θ, b, α) = A*log(max(θ * f(b) - R₂, 0.0))
	P₂(R₂, b, r, s) = s*(max(R̄(b,r) - R₂, 0))
	Obj₂(R₂,b,θ,r,s) = U₂(R₂, θ, b, α) - P₂(R₂, b, r, s)
end

function SPopt(b, r, α, θ, s)
	SP = Model(Ipopt.Optimizer)
	@variable(SP, R₂ ≥ 0)
	@constraint(SP, R₂ ≤ Y₂(b,θ))
	@NLobjective(SP, Max, Obj₂(R₂,b,θ,r,s))
	optimize!(SP)
	return (R₂ = value(R₂), value = objective_value(SP))
end

Analytically, any θ>2.42 should generate R₂ = 120.0

But the following gives R₂ = 93.1506

SPopt(b, r, α, 2.93, s)

and manually, the value function gives a higher value for R₂ = 120

Obj₂(120,b,2.93,r,s)

What am I doing wrong?

odow · February 23, 2024, 7:18pm

Hi @Ayesha, here is what I get when I run your code (I changed @NLobjective to @objective):

julia> using JuMP, Ipopt

julia> begin # set parameters
               b = 100.0
               r = 0.2 
               α = 0.85
               θ = 2.5
               β = 0.9
               A = 1.0 
               s = 1.0 
       end
1.0

julia> begin 
               f(b) = b^α
               Y₂(b,θ) = θ*f(b)
               C₂(R₂,b,θ) = Y₂(b,θ) - R₂
               R̄(b,r) = (1+r)*b
               U₂(R₂, θ, b, α) = A*log(max(θ * f(b) - R₂, 0.0))
               P₂(R₂, b, r, s) = s*(max(R̄(b,r) - R₂, 0))
               Obj₂(R₂,b,θ,r,s) = U₂(R₂, θ, b, α) - P₂(R₂, b, r, s)
       end
Obj₂ (generic function with 1 method)

julia> function SPopt(b, r, α, θ, s)
               SP = Model(Ipopt.Optimizer)
               @variable(SP, 0 <= R₂ <= Y₂(b,θ))
               @objective(SP, Max, Obj₂(R₂,b,θ,r,s))
               optimize!(SP)
               return (R₂ = value(R₂), value = objective_value(SP), summary = solution_summary(SP))
       end
SPopt (generic function with 1 method)

julia> ret = SPopt(b, r, α, 2.93, s);
This is Ipopt version 3.14.14, running with linear solver MUMPS 5.6.2.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        1

Total number of variables............................:        1
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        1
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0 -1.1500067e+02 0.00e+00 9.93e-01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -1.1489258e+02 0.00e+00 5.28e-01  -1.0 1.09e-01    -  5.26e-01 1.00e+00f  1
   2 -1.1320704e+02 0.00e+00 1.69e+00  -1.0 1.70e+00    -  1.00e+00 1.00e+00f  1
   3 -1.1245211e+02 0.00e+00 9.98e-01  -1.0 7.60e-01    -  7.39e-01 1.00e+00f  1
  ... lots of lines omitted ...
2999  3.2901247e+00 0.00e+00 3.73e-02  -2.5 1.79e-01  -0.6 1.00e+00 1.56e-02f  7
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
3000  3.2901843e+00 0.00e+00 3.73e-02  -2.5 8.19e-01  -1.1 1.00e+00 1.95e-03f 10

Number of Iterations....: 3000

                                   (scaled)                 (unscaled)
Objective...............:  -3.2901843244621571e+00    3.2901843244621571e+00
Dual infeasibility......:   3.7325414981530358e-02    3.7325414981530358e-02
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.8477234972151210e-03    2.8477234972151210e-03
Overall NLP error.......:   3.7325414981530358e-02    3.7325414981530358e-02


Number of objective function evaluations             = 55384
Number of objective gradient evaluations             = 3001
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 3000
Total seconds in IPOPT                               = 2.689

EXIT: Maximum Number of Iterations Exceeded.
(R₂ = 120.0000475444581, value = 3.290184324462157, summary = * Solver : Ipopt

* Status
  Result count       : 1
  Termination status : ITERATION_LIMIT
  Message from the solver:
  "Maximum_Iterations_Exceeded"

* Candidate solution (result #1)
  Primal status      : UNKNOWN_RESULT_STATUS
  Dual status        : UNKNOWN_RESULT_STATUS
  Objective value    : 3.29018e+00
  Dual objective value : 1.15175e-02

* Work counters
  Solve time (sec)   : 2.68816e+00
  Barrier iterations : 3000

This finds R₂ = 120, but Ipopt terminated with UNKNOWN_RESULT_STATUS, so anything could have happened. It’s likely something similar happened for you.

The problem is that Ipopt expects problems to be convex and differentiable.

You could get rid of the non-differentiable max with something like:

function SPopt2(b, r, α, θ, s)
    SP = Model(Ipopt.Optimizer)
    @variable(SP, 0 <= R₂ <= Y₂(b,θ))
    @variable(SP, tmp1 >= 0)
    @constraint(SP, tmp1 == θ * b^α - R₂)
    @variable(SP, tmp2 >= 0)
    @constraint(SP, tmp2 >= (1+r)*b - R₂)
    @objective(SP, Max, A*log(tmp1) - s*tmp2)
    optimize!(SP)
    return (R₂ = value(R₂), value = objective_value(SP), summary = solution_summary(SP))
end
ret = SPopt2(b, r, α, 2.93, s)

This yields the much nicer

julia> ret = SPopt2(b, r, α, 2.93, s)
This is Ipopt version 3.14.14, running with linear solver MUMPS 5.6.2.

Number of nonzeros in equality constraint Jacobian...:        2
Number of nonzeros in inequality constraint Jacobian.:        2
Number of nonzeros in Lagrangian Hessian.............:        1

Total number of variables............................:        3
                     variables with only lower bounds:        2
                variables with lower and upper bounds:        1
                     variables with only upper bounds:        0
Total number of equality constraints.................:        1
Total number of inequality constraints...............:        1
        inequality constraints with only lower bounds:        1
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0 -4.6151712e+00 1.47e+02 1.20e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  5.9686771e-01 1.43e+02 1.84e+02  -1.0 7.36e+01    -  1.35e-04 2.52e-02f  1
   2  3.1411827e+00 0.00e+00 2.20e+00  -1.0 1.18e+02    -  1.68e-02 1.00e+00h  1
   3  3.0359454e+00 0.00e+00 1.66e-02  -1.0 1.05e-01    -  9.83e-01 1.00e+00f  1
   4  9.8428081e-01 0.00e+00 7.57e-04  -1.7 2.10e+00    -  9.67e-01 1.00e+00f  1
   5  5.0606236e-02 0.00e+00 2.50e-05  -3.8 9.61e-01    -  9.75e-01 1.00e+00f  1
   6 -2.6511132e-01 0.00e+00 1.08e-06  -5.7 3.25e+00    -  9.61e-01 1.00e+00f  1
   7 -4.7353848e-02 0.00e+00 1.44e-09  -8.6 2.25e-01    -  9.99e-01 1.00e+00f  1
   8  3.1300046e+00 0.00e+00 2.61e-09  -8.6 1.51e+01    -  1.00e+00 2.17e-01h  1
   9  2.5337834e+00 0.00e+00 2.56e-13  -8.6 6.19e-01    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  3.2494826e+00 2.84e-14 1.01e-11 -12.9 3.03e+00    -  9.07e-01 9.87e-01h  1
  11  3.2888312e+00 0.00e+00 6.01e-13 -12.9 1.03e+00    -  1.00e+00 1.00e+00f  1
  12  3.2901615e+00 2.84e-14 6.61e-16 -12.9 3.57e-02    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 12

                                   (scaled)                 (unscaled)
Objective...............:  -3.2901615021128144e-08    3.2901615021128143e+00
Dual infeasibility......:   6.6124070239343904e-16    6.6124070239343909e-08
Constraint violation....:   2.8421709430404007e-14    2.8421709430404007e-14
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   1.2612801725865069e-13    1.2612801725865069e-05
Overall NLP error.......:   1.2612801725865069e-13    1.2612801725865069e-05


Number of objective function evaluations             = 13
Number of objective gradient evaluations             = 13
Number of equality constraint evaluations            = 13
Number of inequality constraint evaluations          = 13
Number of equality constraint Jacobian evaluations   = 1
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations             = 12
Total seconds in IPOPT                               = 0.006

EXIT: Optimal Solution Found.
(R₂ = 120.00031072467335, value = 3.2901615021128143, summary = * Solver : Ipopt

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "Solve_Succeeded"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 3.29016e+00
  Dual objective value : 1.00003e+00

* Work counters
  Solve time (sec)   : 6.28018e-03
  Barrier iterations : 12

Ayesha · February 24, 2024, 11:38pm

Thank you so much for the explanation. The second approach is fast!

What is the reason for using >= instead of == in the constraint of the second approach @constraint(SP, tmp2 >= (1+r)*b - R₂)?

odow · February 25, 2024, 12:34am

I’m using this trick: Tips and tricks · JuMP