# Choice variable doesn't optimize the objective function

**URL:** <https://discourse.julialang.org/t/choice-variable-doesnt-optimize-the-objective-function/110633>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, optim\
**Created:** [February 23, 2024, 10:21am UTC](https://discourse.julialang.org/t/choice-variable-doesnt-optimize-the-objective-function/110633 "2024-02-23T10:21:35Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ayesha](https://avatars.discourse-cdn.com/v4/letter/a/c0e974/32.png) [@Ayesha](https://discourse.julialang.org/u/Ayesha)\
**Post date:** [February 23, 2024, 10:21am UTC](https://discourse.julialang.org/t/choice-variable-doesnt-optimize-the-objective-function/110633/1 "2024-02-23T10:21:35Z")

</div>

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

```

```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

```

```julia
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

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

```

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

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

```

What am I doing wrong?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [February 23, 2024, 7:18pm UTC](https://discourse.julialang.org/t/choice-variable-doesnt-optimize-the-objective-function/110633/2 "2024-02-23T19:18:19Z")

</div>

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

```julia
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:

```julia
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
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

```

---

<div class="post-metadata">

**Author:** ![Ayesha](https://avatars.discourse-cdn.com/v4/letter/a/c0e974/32.png) [@Ayesha](https://discourse.julialang.org/u/Ayesha)\
**Post date:** [February 24, 2024, 11:38pm UTC](https://discourse.julialang.org/t/choice-variable-doesnt-optimize-the-objective-function/110633/3 "2024-02-24T23:38:03Z")

</div>

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₂)`?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [February 25, 2024, 12:34am UTC](https://discourse.julialang.org/t/choice-variable-doesnt-optimize-the-objective-function/110633/4 "2024-02-25T00:34:31Z")

</div>

I’m using this trick: [Tips and tricks · JuMP](https://jump.dev/JuMP.jl/stable/tutorials/linear/tips_and_tricks/#Max)