Problem with cost function in JuMP

Hello everyone,

I am interested in the following optimization problem that I would like to solve with JuMP:

\begin{split} \min_{x \in \mathbb{R}^n} & (z^Tx)^2 - \log(b^T x[l:\mbox{end}])\\ \mbox{ subject to }& x[l:\mbox{end}] \geq 0 \end{split}

with k<n, b \in \mathbb{R}^k and z \in \mathbb{R}^n. x[l:\mbox{end}] denotes the last l components of the vector x.

I am not sure of the way to provide parameters to a user-defined function. Also, I got an error saying that I should use only scalar expressions which seems to be the case

Here is my code:

Code for the cost function

function _cost(x, wquad, w∂,p)
    J = 0.0
    nx = length(x)
    Ne = size(wquad,2)
    for i=1:Ne
        J += 0.5*sum(wquad[j,i]*x[j] for j=1:nx)^2

        if p==0
        J += -log(w∂[1,i]*x[end])
        else
        J += -log(sum(w∂[j,i]*x[end-(p+2)+j] for j=1:p+2))
        end
    end
    J *=  (1/Ne)
    return J
end

Rest of the code for the optimization

p=2
wq = randn(32,200)
w∂ = randn(4,200)

Ne = 200
nx = size(wq,1)
nlog = size(w∂,1)


# Create the model 
model = Model(Ipopt.Optimizer)

# Create variable x
@variable(model, x[1:nx], start = 0.001)

@NLparameter(model, p_order[1]== p)
@NLparameter(model, z[i = 1:nx, j = 1:Ne] == max(wq[i,j], 0.001))
@NLparameter(model, b[i = 1:nlog, j = 1:Ne] == max(w∂[i,j], 0.001))


# Define cost function
cost(x...) = _cost(x, z, b, p_order)
register(model, :cost, 1, cost, autodiff=true)

# Set lower bound on x
for i in nx-nlog:nx
    set_lower_bound(x[i], 0.0)
end


@NLobjective(model, Min, cost(x))

@time optimize!(model)

Code return message

Unexpected array VariableRef[x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10], x[11], x[12], x[13], x[14], x[15], x[16], x[17], x[18], x[19], x[20], x[21], x[22], x[23], x[24], x[25], x[26], x[27], x[28], x[29], x[30], x[31], x[32]] in nonlinear expression. Nonlinear expressions may contain only scalar expressions.

Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] _parse_NL_expr_runtime(::Model, ::Array{VariableRef,1}, ::Array{JuMP._Derivatives.NodeData,1}, ::Int64, ::Array{Float64,1}) at /Users/me/.julia/packages/JuMP/CZ8vV/src/parse_nlp.jl:207
 [3] top-level scope at /Users/me/.julia/packages/JuMP/CZ8vV/src/parse_nlp.jl:82
 [4] top-level scope at /Users/me/.julia/packages/JuMP/CZ8vV/src/macros.jl:1274
 [5] top-level scope at In[37]:28

If you do this, then (I think), x[1] will passed as x, x[2] as z, x[3] as b and x[4:end] as p_order. Is that what you want?

1 Like

Why not just write it out in the macro?

@NLobjective(model, Min, sum(z[i] * x[i] for i=1:N)^2 - log(sum(b[i] * x[i] for i=l:N)))
1 Like

Thanks I think that’s the way to go. I’ve just found easier to write functions called by others functions.