quadratic equality constraint: Julia+Jump+Gurobi

Question

i'm a newbie to julia and struggling how to handle quadratic equality constraints when using gurobi as the main solver. Can u may take a quick look at the following listing? I alread know that such structures are impossible to solve by using gurobi, but how do I can bypass this difficulty and thus enforce the condition to my model?

[listing][1]

Relevant and defined Variables causing problems:

• x: binary decision variable
• s: s>=0

I would be glad about any help!

Andre

eukl_distance = [0 6 5 7 9 8; 6 0 6 3 14 9; 5 6 0 5 9 4; 7 3 5 0 13 8; 
9 14 9 13 0 7; 8 9 4 8 7 0]

# Nachfrage an Knoten i
d = [0,2,1,1,3,1]

no_nodes = 6

# Zeitfenster 
time_windows = [Vector{Float64}(2) for i=1:no_nodes]
time_windows[1] = [0,10]
time_windows[2] = [1,4]
time_windows[3] = [2,5]
time_windows[4] = [4,7]
time_windows[5] = [3,6]
time_windows[6] = [5,8]

# Fahrtzeit von i nach j und Servicezeit an Knoten i
t = [Array{Float64}(no_nodes) for i=1:no_nodes]
t[1] = [0,2,1.5,2.5,4,3.5]
t[2] = [2,0,2,1,6,4]
t[3] = [1.5,2,0,1.5,4,1.25]
t[4] = [2.5,1,1.5,0,5.75,3.5]
t[5] = [4,6,4,5.75,0,2.5]
t[6] = [3.5,4,1.25,3.5,2.5,0]

# Anzahl Fahrzeuge
k = 2

# Kapazität Fahrzeuge
C = 5

# Strafkostenfaktor
z = 2

M = 100000000000000000

using JuMP, Gurobi, LightGraphs, Distances

m = Model(solver = GurobiSolver())

@variable(m, x[i=1:no_nodes,j=1:no_nodes,kk=1:k], Bin)
@variable(m, y[i=1:no_nodes, kk=1:k], Bin)
@variable(m, 0<=s[i=1:no_nodes,kk=1:k]<=M)
@variable(m, w1[i=1:no_nodes, kk=1:k]>=0)
@variable(m, w2[i=1:no_nodes, kk=1:k]>=0)

# Hilfsvariablen
@variable(m, h1[i=1:no_nodes, j=1:no_nodes, kk=1:k]>=0)
@variable(m, h2[i=1:no_nodes, j=1:no_nodes, kk=1:k]>=0)

# Beachte: Division von z durch no_nodes notwendig,
# da ansonsten w1 und w2 für jedes i über jedes j
# addiert werden. Beispiel: Sei i = 2 und j=1:6
# x[2,1,1]*dist+w1[2,1]*2+w2[2,1]*2+x[2,1,2]*dist+w1[2,2]*2+w2[2,2]*2 +
# x[2,2,1]*dist+w1[2,1]*2+w2[2,1]*2+x[2,2,2]*dist+w1[2,2]*2+w2[2,2]*2 + 
# ... x[2,6,1]*dist+w1[2,1]*2+w2[2,1]*2+x[2,6,2]*dist+w1[2,2]*2+w2[2,2]*2
@objective(m, Min, sum(x[i=ii,j=jj,n=kk]*eukl_distance[ii,jj]+w1[i=ii,n=kk]*z/no_nodes+w2[i=ii,n=kk]*z/no_nodes for ii=1:no_nodes, jj=1:no_nodes, kk=1:k))

# erlaube verfrühte sowie verspätete Ankunft
# Abfahrt am Depot zum Zeitpunkt 0
@constraint(m, earlyArrival[ii=1:no_nodes,kk=1:k], time_windows[i=ii][1]-s[i=ii,n=kk]-w1[i=ii,n=kk]<=0)
@constraint(m, lateArrival[ii=1:no_nodes,kk=1:k], time_windows[i=ii][2]-s[i=ii,n=kk]+w2[i=ii,n=kk]>=0)

# Eliminiere Kanten von i zu i
# [i,i]=0
@constraint(m, noSelfConnection[ii=1:no_nodes,kk=1:k], x[i=ii,j=ii,n=kk]==0)

# Jeder Kunde wird von einem Fahrzeug besucht
@constraint(m, oneCustomerVisitation[ii=2:no_nodes], sum(y[i=ii,kk=1:k])==1)

# Genau K Fahrzeuge erreichen das Depot
@constraint(m, sum(y[i=1,kk=1:k])==k)

# Wenn ein Fahrzeug einen Kunden i anfährt,
# muss die Summe aller genutzten ausgehenden und eingehenden Kanten,
# basienrend auf Knoten i, genau 1 entsprechen
@constraint(m, completeRoute1[ii=1:no_nodes, kk=1:k], sum(x[i=ii,j=1:no_nodes,n=kk])-y[i=ii,n=kk]==0)
@constraint(m, completeRoute2[ii=1:no_nodes, kk=1:k], sum(x[j=1:no_nodes,i=ii,n=kk])-y[i=ii,n=kk]==0)

# Kapazitätsbeschränkung
@constraint(m, capacityConstr[kk=1:k], sum(d[i]*y[i,kk] for i in 1:no_nodes)<=C)

# Ankunft am Knoten i entspricht der 
# Summe aus Ankunftszeit am Knoten j und Lieferzeit 
# von Knoten i zu j
for i=1:no_nodes, j=1:no_nodes, kk=1:k
    @constraint(m, h1[i,j,kk]<=M*x[i,j,kk])
    @constraint(m, h1[i,j,kk]<=x[i,j,kk]*s[i,kk])
    @constraint(m, h1[i,j,kk]>=s[i,kk]-M*(1-x[i,j,kk]))

    @constraint(m, h2[i,j,kk]<=M*x[i,j,kk])
    @constraint(m, h2[i,j,kk]<=x[i,j,kk]*s[j,kk])
    @constraint(m, h2[i,j,kk]>=s[j,kk]-M*(1-x[i,j,kk]))
    if i==1
        @constraint(m, x[i,j,kk]*t[i][j]-h2[i,j,kk]==0)
    else
        @constraint(m, h1[i,j,kk]+x[i,j,kk]*t[i][j]-h2[i,j,kk]==0)
    end
end

# ------------ Lazy Constraint -------------

function buildGraph(n)
g = DiGraph(n)
for i in 1:n, j in 1:n, kk in 1:k
    if round.(Int,getvalue(x)[i,j,kk])==1
        add_edge!(g,(i,j))
    end
end
return simplecycles_hadwick_james(g)
end

function subtour(cb)
subtours = buildGraph(no_nodes)
for n in 1:length(subtours), kk in 1:k
    if subtours[n][1]!=1 
        @lazyconstraint(cb, sum(x[i=subtours[n],j=subtours[n],kk])<=length(subtours[n])-1)
    end
end
end

addlazycallback(m, subtour)

solve(m)
getobjectivevalue(m)

for i in 1:no_nodes, j in 1:no_nodes, kk in 1:k
    if round.(Int,getvalue(x)[i,j,kk])==1
        println("($i,$j,$kk) = 1")
    end
end