Periodic cost function for JuMP

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.

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:

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.

You cannot use these sorts of objective functions directly in JuMP. Read:

Instead, you need to formulate it as a mixed-integer quadratic program.

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)

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.

Okay, I tried your suggested code with Gurobi and get the following output:

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:
image

Are there any options for Gurobi that I should tune for this type of problem?
Thanks very much!

I probably made a mistake in the formulation. I’ll take another look.

Oops. The model I originally wrote was for sum(@.(max(..., not sum(@.(min(....

Perhaps this will work

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)
2 Likes

Thanks, this works perfectly!

1 Like