# Unexpected array VariableRef[c c c] in nonlinear expression

Code:

``````function b̂ₜfunc(
x̂ₜ₊₁::AbstractVector,
Σ̂ₜ₊₁::AbstractMatrix,
γ::AbstractFloat,
ξ::AbstractFloat,
b̂ₜ::AbstractVector
)
n_assets = length(x̂ₜ₊₁)
model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
@variable(model, c[1:n_assets])
@constraint(model, sum(c)==0)
@NLobjective(model, Min, γ*c'*Σ̂ₜ₊₁*(c')' + c'*((2γ*b̂ₜ'*Σ̂ₜ₊₁)'-x̂ₜ₊₁) + ξ*sum(abs(item) for item in c))
optimize!(model)
return value.(c)
end
``````

Reproducible example:

``````x̂ₜ₊₁ = [0.9884, 1.02564, 1.01561];
Σ̂ₜ₊₁ = rand(3, 3);
gam = 0.02;
kasi = 0.03;
b = [0.5, 0.1, 0.4];
b̂ₜfunc(x̂ₜ₊₁, Σ̂ₜ₊₁, gam, kasi, b)
``````

Throws:

``````ERROR: Unexpected array VariableRef[c c c] in nonlinear expression. Nonlinear expressions may contain only scalar expressions.
``````

Other variants that I’ve tried:

``````@NLobjective(model, Min, sum((γ*c'*Σ̂ₜ₊₁)[item]*c[item] for item=1:n_assets) + c'*((2γ*b̂ₜ'*Σ̂ₜ₊₁)'-x̂ₜ₊₁) + ξ*sum(abs(item) for item in c))

# And

@NLobjective(model, Min, γ*c'*Σ̂ₜ₊₁*(c')' + sum(c[item]*((2γ*b̂ₜ'*Σ̂ₜ₊₁)'-x̂ₜ₊₁)[item] for item=1:n_assets) + ξ*sum(abs(item) for item in c))

# And
@NLobjective(model, Min, γ*c'*Σ̂ₜ₊₁*c + only(c'*((2γ*b̂ₜ'*Σ̂ₜ₊₁)'-x̂ₜ₊₁)) + ξ*sum(abs(item) for item in c))
``````

Throws the same error.

When I try `@NLobjective(model, Min, γ*c'*Σ̂ₜ₊₁*c)` I get the same error. Hence, the problem is with the matrix multiplication. I do not understand. This is the same as this example. Why it does not work in my case? Got it. Because I used `@NLobjective`.

The problem was with `@NLobjective`. I replaced it with `@objective`. Also, there is no need to use `only`. Furthermore, `sum(abs(item) for item in c)` can be replaced with `sum(abs.(c))`. Final functions:

``````function b̂ₜfunc(
x̂ₜ₊₁::AbstractVector,
Σ̂ₜ₊₁::AbstractMatrix,
γ::AbstractFloat,
ξ::AbstractFloat,
b̂ₜ::AbstractVector
)
n_assets = length(x̂ₜ₊₁)
model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
@variable(model, c[1:n_assets])
@constraint(model, sum(c)==0)
@objective(model, Min, γ*c'*Σ̂ₜ₊₁*c + c'*((2γ*b̂ₜ'*Σ̂ₜ₊₁)'-x̂ₜ₊₁) + ξ*sum(abs.(c)))
optimize!(model)
return value.(c)
end

x̂ₜ₊₁ = [0.9884, 1.02564, 1.01561];
Σ̂ₜ₊₁ = rand(3, 3);
gam = 0.02;
kasi = 0.03;
b = [0.5, 0.1, 0.4];
b̂ₜfunc(x̂ₜ₊₁, Σ̂ₜ₊₁, gam, kasi, b)
# 3-element Vector{Float64}:
#   1.4843778267918864e20
#  -6.863558855303889e18
#  -1.4157422382388475e20
``````

Wow Such big values.

1 Like

You should always check the `termination_status(model)` before querying values. Did the problem solve correctly?

And yes, the `@objective` macro now supports nonlinear and it fixes a number of limitations from the `@NL` macros.

1 Like

If this can be implied by the function that you mentioned:

``````julia> termination_status(model)
NORM_LIMIT::TerminationStatusCode = 17
``````

It’s not clear. I do not know what’s the meaning of `NORM_LIMIT::TerminationStatusCode = 17`. I guess it’s not converged. I don’t know.

Ah, Grrreat!! Thank you so much. `NORM_LIMIT` means that the solver did not terminate with an optimal solution. See JuMP · JuMP

If you turn on printing, you’ll see Ipopt print:

`EXIT: Iterates diverging; problem might be unbounded.`

The issue is that Ipopt assumes problems are smooth and twice differentiable. `abs(c)` violates that assumption, so Ipopt is not guaranteed to find a solution.

One option is the replace the L1 penalty `sum(abs.(c))` with the L2 `sum(c.^2)`.

Another is to use a reformulation (see documentation Tips and tricks · JuMP):

``````using JuMP, Ipopt
function b̂ₜfunc(
x̂ₜ₊₁::AbstractVector,
Σ̂ₜ₊₁::AbstractMatrix,
γ::AbstractFloat,
ξ::AbstractFloat,
b̂ₜ::AbstractVector
)
n_assets = length(x̂ₜ₊₁)
model = Model(Ipopt.Optimizer)
@variable(model, c[1:n_assets])
@constraint(model, sum(c) == 0)
@variable(model, t[1:n_assets] >= 0)
@constraint(model, t .>= c)
@constraint(model, t .>= -c)
@objective(model, Min, γ*c'*Σ̂ₜ₊₁*c + c'*((2γ*b̂ₜ'*Σ̂ₜ₊₁)'-x̂ₜ₊₁) + ξ*sum(t))
optimize!(model)
return value.(c)
end

x̂ₜ₊₁ = [0.9884, 1.02564, 1.01561];
Σ̂ₜ₊₁ = rand(3, 3);
γ = 0.02;
ξ = 0.03;
b̂ₜ = [0.5, 0.1, 0.4];
b̂ₜfunc(x̂ₜ₊₁, Σ̂ₜ₊₁, γ, ξ, b̂ₜ)
``````
1 Like

Thank you so much. I shall rely much more on the provided documentation. Thanks.

Since I’m trying to implement a proposed algorithm, I should follow the same signature (Although I’m aware that it does not change the result of the model). However, I got your point on this. Thank you so much. This is the most interesting part! However, it’s a little bit ambiguous for me. For example, you’ve set `t[1:n_assets] >= 0`, and later you noted `t .>= -c`. Isn’t this redundant? Or I’m missing something here. Oh, I got it. Because we do not know whether c_i would be a positive or negative value. Hence, we need both `t .>= c` and `t .>= -c` to ensure that t_i would result in the absolute value of c_i. Thank you so much! It’s a beautiful trick.

1 Like

Yes, technically the `t >= 0` is redundant. But it doesn’t hurt to have it.

1 Like