Say you have *n* different shops in which you want to launch a product for *n* different months, and you want to maximize profits. The problem is determining the order of the launch sequence of the product when the price in shop_*i* refers to the average price for the product in some given set of the shops. It is only possible to launch in one shop for every month, and it is only possible to launch one time in each shop. For every shop there is an initial price of the product, some given customer size, and baskets referring prices. Look at the following example:

initial_price = [1,2,4],

customer_size = [20, 15, 10]

referring_basket_shop_1 = {price_shop_2,price_shop_3}

referring_basket_shop_2 = {price_shop_3}

referring_basket_shop_3 = {price_shop_1, price_shop_2}

The first shop to launch the product initiates the price with the initial_price parameter. If referring_basket_shop_*i* contains no launched shops, price_shop_*i* = initial_price_*i*. If referring_basket_shop_*i* contains launched shops, price_shop_*i* will be the average price of the launched prices in referring_basket_shop_*i*.

Calculating all possible outcomes and determining the best one is a *n!* problem, so I want to use Julia optimization (I donâ€™t know which solver is the best, but in the following example I used Ipopt).

I **tried** to formulate the problem like this:

**x** is a *n* by *n* matrix containing the price for shop_i and month_i for i, j = 1:*n*.

**y** is a *n* by *n* matrix containing binary elements as decision variables i.e. if shop_*i* launches in month_j for i,j = 1:*n*.

**nr** is a *n* by *n* matrix containing only integer elements and itâ€™s used to calculate the average i.e. ensuring that we divide with the correct number.

Here is my code:

```
using JuMP
using Ipopt
using MathOptInterface
using Juniper
n = 3
optimizer = Juniper.Optimizer
nl_solver = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
m = Model(optimizer_with_attributes(optimizer, "nl_solver"=>nl_solver))
@variable(m, x[i = 1:n, j = 1:n] >= 0)
@variable(m, y[i = 1:n, j = 1:n], Bin)
@variable(m, nr[i = 1:n, j = 1:n], Int)
start = [1,2,4]
pop = [20,15,10]
@objective(m, Max, sum(x[i,j]*y[i,j] * pop[j] for i =1:n, j = 1:n) )
#Launch one time per shop
for j = 1:n
@constraint(m, sum(y[i,j] for i = 1:n) == 1 )
end
#Launch one time per month
for i = 1:n
@constraint(m, sum(y[i,j] for j = 1:n) == 1 )
end
#Initial price
for i = 1:n
@constraint(m, x[i,1]*y[i,1] <= start[i] )
end
#Constraining nr
for i = 1:n
for j = 1:n
@constraint(m, nr[i,j] >= 1 )
@constraint(m, nr[i,j] <= n )
end
end
#Calculating launch sequence
for j = 2:n
@constraint(m, nr[1,j] >= sum(y[2,i] + y[3,i] for i = 1:j-1) )
@NLconstraint(m, x[1,j] == sum(y[2,i]*x[2,i] + y[3,i]*x[3,i] for i = 1:j-1) / nr[1,j] )
@constraint(m, nr[2,j] >= sum(y[3,i] for i = 1:j-1) )
@NLconstraint(m, x[2,j] == sum(y[3,i]*x[3,i] for i = 1:j-1) / nr[2,j] )
@constraint(m, nr[3,j] >= sum(y[2,i] + y[1,i] for i = 1:j-1) )
@NLconstraint(m, x[3,j] == sum(y[2,i]*x[2,i] + y[1,i]*x[1,i] for i = 1:j-1) / nr[3,j] )
end
optimize!(m)
println(JuMP.value.(x))
println(JuMP.value.(y))
println(JuMP.value.(nr))
println(JuMP.objective_value(m))
println(JuMP.termination_status(m))
```

My code does not give any useful results, and I am stuck. Any suggestions?