# Interpolation in constraint

**URL:** <https://discourse.julialang.org/t/interpolation-in-constraint/7325>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [November 26, 2017, 6:28pm UTC](https://discourse.julialang.org/t/interpolation-in-constraint/7325 "2017-11-26T18:28:55Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![MaximilianJHuber](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maximilianjhuber/32/2579_2.png) [@MaximilianJHuber](https://discourse.julialang.org/u/MaximilianJHuber)\
**Post date:** [November 26, 2017, 6:28pm UTC](https://discourse.julialang.org/t/interpolation-in-constraint/7325/1 "2017-11-26T18:28:55Z")

</div>

Hi, I have a conceptional question:

I have 10 non-linear constraints for k=k\_1, ..., k\_{10} (that represent a discrete grid of a continuous variable) and 10 coefficients, \theta = \theta\_1, ... \theta\_{10}, that are determined by those:

u'\Big(\Psi(k) \cdot \theta, \gamma\Big) =\beta \, u'\Big( \Psi \big( f(k, \alpha, A) +(1- \delta)k- \Psi(k) \cdot \theta\big)\cdot \theta , \gamma\Big) \cdot \quad \quad \quad\quad\quad\quad\quad\quad\quad\Big( f'\big(f(k, \alpha, A) +(1- \delta)k - \Psi(k) \cdot \theta, \alpha, A\big) + (1-\delta) \Big)

\Psi evaluates 10 basis functions (in my case of a Chebyshev polynomials) used to approximate a univariate function.  
While \Psi(k) can be evaluated before assembling the JuMP model, \Psi \big( f(k, \alpha, A) +(1- \delta)k- \Psi(k) \cdot \theta\big) is a function of the variables of the JuMP model: \beta, \alpha, A, \gamma, \delta, \theta[1:10].

I see two possible approaches:

- Register \Psi in JuMP and hope AutoDiff can handle it:  
Chebyshev polynomials are [defined](https://en.wikipedia.org/wiki/Chebyshev_polynomials#Definition) iteratively, I believe AutoDiff should be able to do this in theory. But, using the [BasisMatrices](https://github.com/QuantEcon/BasisMatrices.jl) package throws an error:

```julia
MethodError: no method matching evalbasex!(::Array{Any,2}, ::Array{JuMP.NonlinearExpression,1}, ::BasisMatrices.ChebParams{Float64}, ::Array{JuMP.NonlinearExpression,1})
Closest candidates are:
  evalbasex!(::AbstractArray{T,2} where T, ::AbstractArray{T<:Number,1}, ::BasisMatrices.ChebParams, ::AbstractArray{T<:Number,1}) where T<:Number at /home/juser/.julia/v0.6/BasisMatrices/src/cheb.jl:168
  evalbasex!(::AbstractArray{T,2} where T, ::BasisMatrices.ChebParams, ::Any) at /home/juser/.julia/v0.6/BasisMatrices/src/cheb.jl:194

```

- Update the constraints every iteration of the optimization:  
It is easy to calculate f(k, \alpha, A) +(1- \delta)k- \Psi(k) \cdot \theta every iteration, if I could extract it from the JuMP model and evaluate \Psi there, I could update the constraint. This is discussed [here](https://discourse.julialang.org/t/updating-data-in-jump-constraints-iteratively/5467) and [here](http://www.juliaopt.org/JuMP.jl/0.18/probmod.html#modifying-constraints) and not yet supported by JuMP.

Should I code up \Psi as “pure” as possible myself, use a different package, or wait for the updating functionality?

Here is the full model:

```julia
k_stst = 4.628988089138438 #reference point of the grid
basis = Basis(ChebParams(10, 0.2*k_stst, 2*k_stst))
K = nodes(basis)[1] #grid
Ψ = BasisMatrix(basis, Expanded(), K).vals[1] #\Psi(k)

function f(k, α, A)
    A*k^α
end

function f_prime(k, α, A)
    A*α*k^(α-1)
end

u_crra_prime(c, γ) = c.^-γ

```

```julia
@variable(m, β, start = 1)
@variable(m, δ, start = 1)
@variable(m, α, start = 1)
@variable(m, A, start = 1)
@variable(m, γ, start = 1)

@variable(m, θ[1:10], start = 0)
setvalue(θ[1], 0.001)

JuMP.register(m, :u_crra_prime, 2, u_crra_prime, autodiff=true)
JuMP.register(m, :f, 3, f, autodiff=true)
JuMP.register(m, :f_prime, 3, f_prime, autodiff=true)

@NLexpression(m, Kprime[i=1:10], f(K[i], α, A) + (1-δ)*K[i] - sum(Ψ[i, k] * θ[k] for k in 1:10)) #argument of \Psi

```

```julia
function Ψprime(k)
    return BasisMatrix(Basis(ChebParams(10, 0.2*4.628988089138438, 2*4.628988089138438)), Expanded(), [k]).vals[1]
end
JuMP.register(m, :Ψprime, 1, Ψprime, autodiff=true)

@NLconstraint(m, EE[i=1:10], u_crra_prime(sum(Ψ[i, k] * θ[k] for k in 1:10), γ) == 
    β*u_crra_prime(sum(Ψprime(Kprime[i])[k] * θ[k] for k in 1:10), γ) * 
    (f_prime(Kprime[i], α, A)))

```

The objective is a GMM-type function that compares predicted values against data:

```julia
@NLobjective(m, Min, sum((dataC[t] - predictedC[t])^2 for t in 1:100) + 
    sum((dataK[t] - predictedK[t])^2 for t in 2:100))

```

---

<div class="post-metadata">

**Author:** ![sglyon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sglyon/32/67_2.png) [@sglyon](https://discourse.julialang.org/u/sglyon)\
**Post date:** [November 28, 2017, 5:04pm UTC](https://discourse.julialang.org/t/interpolation-in-constraint/7325/2 "2017-11-28T17:04:00Z")

</div>

Looks like this is a problem with BasisMatrices. Will you please open an issue there describing the issue?

---

<div class="post-metadata">

**Author:** ![MaximilianJHuber](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maximilianjhuber/32/2579_2.png) [@MaximilianJHuber](https://discourse.julialang.org/u/MaximilianJHuber)\
**Post date:** [November 28, 2017, 5:55pm UTC](https://discourse.julialang.org/t/interpolation-in-constraint/7325/3 "2017-11-28T17:55:50Z")

</div>

I opened an [issue](https://github.com/QuantEcon/BasisMatrices.jl/issues/44).  
Still I would be grateful for an educated guess or some advice on evaluating a basis function, that has a for loop in it, inside of an constraint.

---

<div class="post-metadata">

**Author:** ![sglyon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sglyon/32/67_2.png) [@sglyon](https://discourse.julialang.org/u/sglyon)\
**Post date:** [December 7, 2017, 3:59pm UTC](https://discourse.julialang.org/t/interpolation-in-constraint/7325/4 "2017-12-07T15:59:17Z")

</div>

Great, thank you.

I have posted a potential solution to problem on that issue. Could you please give it a try and let me know how it goes?

I only showed code for a working example using ForwardDiff.jl, but would love to see one that uses Jump!