Using JuMP and Ipopt I cannot get a simple Euler constrained optimisation problem to output correct answers

Specific Domains Optimization (Mathematical)
1

Here is the code

using JuMP
using Ipopt

function Euler(x...)
  beta = 0.9
  t = 2
  x = [i for i in x]
  logx = log.(x)
  answer = -(beta .^ ((transpose(1:1:t) .- 1) * logx))
  return answer[1] 
end

m = Model(Ipopt.Optimizer)
@variable(m, x[1:2] >= eps())
@NLconstraint(m, x[1]+x[2] <=20)
register(m, :Euler, 2, Euler, autodiff=true)
optimize!(m)
@show value.(x)

In a two period Euler equation model with an endowment of 20 and a discount rate of 0.9 the maximum utility obtainable is with a consumption bundle of roughly [10.5,9.5]. Despite this I get some weird interior point where consumption is the same in both periods [2.5,2.5]. I tried constraining the problem so that the inequality constraint becomes an equality constraint but then Julia gives me a bundle of [10,10]. I am not sure what I am doing wrong

1 Like
2

Are you missing some code? You don’t seem to actually use the registered function.

3

I edited my original code snippet because I forgot the statements about importing the relevant Julia libraries. My user created function is called “Euler” and that is the function that I use the JuMP register function for. The code executes on my machine so I am confused about there being something missing

1 Like
4

As was pointed out by odow, although you define and register your Euler function, it is not used anywhere in the optimization problem itself.
In particular, you do not define any objective: the problem that Ipopt sees is

\begin{array}{rl} \displaystyle \min_{x_1, x_2} \quad & 0 \\ s.t. \quad & x_1 + x_2 \leq 20,\\ & x_{1}, x_{2} \geq \epsilon. \end{array}

Therefore, Ipopt may return any solution (x_1, x_2) such that x_1, x_2 > \epsilon and x_1 + x_2 \leq 20.

If you want to minimize Euler(x), you need to declare an objective with @NLObjective, see e.g. this example

2 Likes
5

Hello there thank you for the reply. What happened is that I transcribed the code from someone else and that error I carelessly carried over.

My issue persists even when I set out the line of code that runs the internal solver. I can get an answer but it is an answer which doesn’t make sense from the perspective of an Euler equation where future consumption is valued less than present consumption. The inequality constraint is not used optimally because consumption in both periods doesn’t equal the budget.

using JuMP
using Ipopt

function Euler(x...)
  beta = 0.9
  t = 2
  x = [i for i in x]
  logx = log.(x)
  answer = -(beta .^ ((transpose(1:1:t) .- 1) * logx))
  return answer[1] 
end

##

m = Model(Ipopt.Optimizer)
@variable(m, x[1:2] >= eps())
@NLconstraint(m, x[1]+x[2] <=20)
register(m, :Euler, 2, Euler, autodiff=true)
@NLobjective(m, Min, Euler(x...))
# set_start_value(x[1], 10.0) setting starting values gives me errors
# set_start_value(x[2], 10.0) setting starting values gives me errors
optimize!(m)
@show value.(x)

I don’t know if you are familiar with MATLAB but what is inexplicable to me is why I can get the solver to converge perfectly on the right answer every time but I am having so much difficulty using Julia for the same basic problem

function V = FlowUtility(T, Beta, C)
%
% PURPOSE: calculates the total utility of consumption assuming an
% additively separable utility function and discount rate Beta
%
% INPUTS: C : Tx1 vector of independent variable
%         T : scalar time
%         Beta : scalar discount rate
%
% OUTPUT: V : -utility

% c = C(:,1);
t = 1:1:T;
V = Beta.^(t-1)*log(C);
V = -V;

return


Beta = 0.9;
T = 2;
K1 = 20;
lb = eps*ones(2,1);
ub = 20*ones(2,1);
guess = 10*ones(2,1);
A = ones(1,2)
opt = optimset('TolFun', 1E-20, 'TolX',...
    1E-20, 'algorithm','sqp');
c = fmincon(@(C) FlowUtility(T,Beta,C),...
    guess,A,K1,[],[],lb,...
    ub,[],opt)
6

There are a few things going on, the main one being that your Euler function is not the same as Matlab’s. You have the equivalent to Beta.^((t-1).*log(C)), not Beta.^(t-1)*log(C).

Here’s how I would write your model:

using JuMP, Ipopt
beta, t = 0.9, 2
model = Model(Ipopt.Optimizer)
@variable(model, 0.001 <= x[1:t] <= 20, start = 10)
@constraint(model, sum(x) <= 20)
@NLobjective(model, Min, -sum(beta^(i - 1) * log(x[i]) for i in 1:t))
optimize!(model)
7

Thank you Oscar. I didn’t reply initially but your solution helped me to get to grips with how Julia works. I would like to use this function with a non-linear constraint like I would with MATLAB but I am struggling with how to use @NLconstraint compared to how I usually do things with a separate function. Here is the function I have so far

using Ipopt, JuMP, Plots
function dynamicChoice(k, t, beta)

    ## Model estimation

    model = Model(Ipopt.Optimizer)
    @variable(model, 0.001 <= x[1:t] <= k, start = k/t)
    @constraint(model, sum(x) <= k)
    @NLobjective(model, Min, -sum(beta^(i - 1) * log(x[i]) for i in 1:t))
    optimize!(model)
    # @show value.(x)

    ## Plotting

    kvector = k*ones(t,1)
    c = accumulate(+, value.(x))  
    stock =  kvector - c

    yvalues = [stock value.(x)]
    xvalues = 1:t
    # theme(::ggplot2)
    plot(xvalues,
            yvalues,
            xlim=(1,t),
            xticks=(1:1:t),
            ylim=(0,k),
            xlabel = "Period",
            ylabel = "Consumption and Investment",
            label = ["capital" "consumption"])
end

I want to have an inequality constraint that includes this for loop

for i in 1:t
   ceq[t] = k[t+1]-Theta*(k[t]-c[t])^Alpha 
end

for given parameters with function defined as

function eulerConstraint(c,t,k,Alpha,Theta)
   return ceq
end

But I am not making much progress

8

I’m not sure I understand the question. Your for loop doesn’t include i, and I don’t know where you use the eulerConstraint function.

Do you want something like:

@NLconstraint(model, [t=1:2], k[t+1]-Theta*(k[t]-c[t])^Alpha <= 0)