# Periodic cost function for JuMP
**URL:** https://discourse.julialang.org/t/periodic-cost-function-for-jump/70760
**Category:** Optimization (Mathematical)
**Created:** [November 1, 2021, 9:11pm UTC](https://discourse.julialang.org/t/periodic-cost-function-for-jump/70760 "2021-11-01T21:11:35Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@sefffal](https://discourse.julialang.org/u/sefffal)
#### Post date: [November 1, 2021, 9:11pm UTC](https://discourse.julialang.org/t/periodic-cost-function-for-jump/70760/1 "2021-11-01T21:11:35Z")
Hello, I am trying to model a tricky optimization problem using JuMP.
**Exposition:**
I am trying to manufacture something a bit like a [Fresnel lens](https://en.wikipedia.org/wiki/Fresnel_lens).
I have a continuous function of height vs x & y position, but we can only manufacture it at 16 discrete levels, and the height of our function is much greater than the highest level we can manufacture.
Thankfully, we are illuminating the surface with monochromatic light of exactly one wavelength. This means, we can “wrap” the surface up and down in increments of one wavelength and it will work identically.
**What I tried:**
My first approach was to create a matrix variable in jump, and create a cost function using `mod`:
```julia
model = Model(COSMO.Optimizer)
@variable(model, z[1:21, 1:21], lower_bound=0, upper_bound=16)
periodic_euclidean_cost(z,z′) = sum(@.(
min(
mod((z - z′)^2, 16),
16 - mod((z - z′)^2, 16)
)^2))
@objective(model, Min, periodic_euclidean_cost(z, z′))
```
My goal was to treat each variable as periodic, and find the closest value modulo 16.
**The Error I Encountered**
Attempting to add that objective to the model leads to `No method matching mod(::QuadExpr, ::Int64)`.
I guess this is not surprising that JuMP can’t optimize through `mod`.
Is there another way I can specify this problem that JuMP will understand?
Thanks!
N.B. I know it would be easy to just wrap our matrix ourselves and not use JuMP, but I have some additional constraints I need JuMP to optimize e.g. the height of the layers.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [November 1, 2021, 9:49pm UTC](https://discourse.julialang.org/t/periodic-cost-function-for-jump/70760/2 "2021-11-01T21:49:11Z")
You cannot use these sorts of objective functions directly in JuMP. Read:
- [https://jump.dev/JuMP.jl/stable/background/should\_i\_use/#Black-box,-derivative-free,-or-unconstrained-optimization](https://jump.dev/JuMP.jl/stable/background/should_i_use/#Black-box,-derivative-free,-or-unconstrained-optimization)
Instead, you need to formulate it as a mixed-integer quadratic program.
```nohighlight
using JuMP, Gurobi
model = Model(Gurobi.Optimizer)
set_optimizer_attribute(model, "NonConvex", 2)
z′ = rand(0:16, 21, 21)
@variables(model, begin
0 <= z[1:21, 1:21] <= 16, Int
t[1:21, 1:21]
div[1:21, 1:21], Int
0 <= rem[1:21, 1:21] <= 15, Int
end)
@constraints(model, begin
[i=1:21, j=1:21], (z[i, j] - z′[i, j])^2 == 16 * div[i, j] + rem[i, j]
[i=1:21, j=1:21], t[i, j] >= rem[i, j]
[i=1:21, j=1:21], t[i, j] >= 16 - rem[i, j]
end)
@objective(model, Min, sum(t[i, j]^2 for i in 1:21, j in 1:21))
optimize!(model)
```
---
### Author:  [@sefffal](https://discourse.julialang.org/u/sefffal)
#### Post date: [November 1, 2021, 10:43pm UTC](https://discourse.julialang.org/t/periodic-cost-function-for-jump/70760/3 "2021-11-01T22:43:42Z")
Thanks for your response @odow!
This is quite a clever way of specifying the problem. I am applying for an academic Gurobi license and will give it a try.
---
### Author:  [@sefffal](https://discourse.julialang.org/u/sefffal)
#### Post date: [November 1, 2021, 11:18pm UTC](https://discourse.julialang.org/t/periodic-cost-function-for-jump/70760/4 "2021-11-01T23:18:40Z")
Okay, I tried your suggested code with Gurobi and get the following output:
```julia
cademic license - for non-commercial use only - expires 2022-10-27
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 242 rows, 484 columns and 484 nonzeros
Model fingerprint: 0x2b837d1e
Model has 121 quadratic objective terms
Model has 121 quadratic constraints
Variable types: 121 continuous, 363 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
QMatrix range [1e+00, 1e+00]
QLMatrix range [1e+00, 8e+01]
Objective range [0e+00, 0e+00]
QObjective range [2e+00, 2e+00]
Bounds range [2e+01, 2e+01]
RHS range [2e+01, 2e+01]
QRHS range [3e+01, 2e+03]
Presolve removed 138 rows and 276 columns
Presolve time: 0.30s
Presolved: 260 rows, 208 columns, 676 nonzeros
Presolved model has 52 quadratic objective terms
Presolved model has 52 bilinear constraint(s)
Variable types: 52 continuous, 156 integer (0 binary)
Root relaxation: objective 8.917000e+03, 451 iterations, 0.00 seconds
H 3676 3984 11745.000000 8917.00000 24.1% 2.8 1s
H 4017 3019 10557.000000 8917.00000 15.5% 2.8 1s
* 4024 2851 798 10494.000000 8917.00000 15.0% 2.8 1s
H 4601 2389 9864.0000000 8917.00000 9.60% 2.7 1s
H 4601 2269 9801.0000000 9721.00000 0.82% 2.7 1s
Cutting planes:
MIR: 1
Explored 4641 nodes (13520 simplex iterations) in 1.71 seconds
Thread count was 8 (of 8 available processors)
Solution count 6: 9801 9864 10494 ... 11808
Optimal solution found (tolerance 1.00e-04)
Best objective 9.801000000000e+03, best bound 9.801000000000e+03, gap 0.0000%
User-callback calls 9456, time in user-callback 0.00 sec
```
So far so good. But when I look at the output for a simple test case, it does not appear well-converged.
Left is `z′`, middle is manually wrapping the values to within 0:16, and right is `z` after optimization:
Are there any options for Gurobi that I should tune for this type of problem?
Thanks very much!
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [November 2, 2021, 12:17am UTC](https://discourse.julialang.org/t/periodic-cost-function-for-jump/70760/5 "2021-11-02T00:17:15Z")
I probably made a mistake in the formulation. I’ll take another look.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [November 2, 2021, 12:29am UTC](https://discourse.julialang.org/t/periodic-cost-function-for-jump/70760/6 "2021-11-02T00:29:59Z")
Oops. The model I originally wrote was for `sum(@.(max(...`, not `sum(@.(min(...`.
Perhaps this will work
```julia
using JuMP, Gurobi
model = Model(Gurobi.Optimizer)
set_optimizer_attribute(model, "NonConvex", 2)
z′ = rand(0:16, 21, 21)
@variables(model, begin
0 <= z[1:21, 1:21] <= 16, Int
div[1:21, 1:21], Int
-8 <= rem[1:21, 1:21] <= 8, Int
end)
@constraints(model, begin
[i=1:21, j=1:21], (z[i, j] - z′[i, j])^2 == 16 * div[i, j] + rem[i, j]
end)
@objective(model, Min, sum(rem[i, j]^2 for i in 1:21, j in 1:21))
optimize!(model)
```
---
### Author:  [@sefffal](https://discourse.julialang.org/u/sefffal)
#### Post date: [November 2, 2021, 1:57am UTC](https://discourse.julialang.org/t/periodic-cost-function-for-jump/70760/7 "2021-11-02T01:57:22Z")
Thanks, this works perfectly!