We store the network in a structure, where cost is an adjacency matrix storing the cost for each edge in the network, we also store the quality parameters (or requirements) in q and the capacities in cap,
In [1]:
immutable PoolingNetwork
ns :: Int
np :: Int
nt :: Int
cost :: SparseMatrixCSC{Float64,Int}
q :: Vector{Float64}
cap :: Vector{Float64}
end
The following block of code is not important, it just plots a pooling network.
In [2]:
using Compose
immutable Node
x :: Float64
y :: Float64
radius :: Float64
label :: String
text :: String
end
draw_node(n::Node) = compose(context(), ( compose(context(), circle(n.x, n.y, n.radius)),
compose(context(), text(n.x, n.y, n.label, hcenter, vcenter)),
compose(context(), text(n.x, n.y - (2*n.radius), n.text, hcenter, vbottom)) ))
function draw_edge(a::Node, b::Node, label::String)
p = [b.x-a.x, b.y-a.y] ./ norm([b.x-a.x, b.y-a.y])
compose(context(),
(compose(context(), line([(a.x+p[1]*a.radius, a.y+p[2]*a.radius), (b.x-p[1]*b.radius, b.y-p[2]*b.radius)])),
compose(context(), text(a.x+(b.x-a.x)/2, a.y+(b.y-a.y)/2+0.02, label, hcenter, vtop))) )
end
function draw_network(G::PoolingNetwork)
dx, dy = 0.4, 1.0/G.ns
N = [ [ Node(0.1, dy*(i-0.5),.05,string(i),string("q<sub>",i,"</sub>")) for i=1:G.ns ];
[ Node(0.1+dx, dy*(i-0.5),.05,string(i+G.ns),string("w<sub>",i+G.ns,"</sub>")) for i=1:G.np ];
[ Node(0.1+2*dx,dy*(i-0.5),.05,string(i+G.ns+G.np),string("q<sub>",i+G.ns+G.np,"</sub>")) for i=1:G.nt ] ]
p = [draw_node(ni) for ni=N]
for j=1:G.ns+G.np+G.nt
for i=find(G.cost[:,j])
push!(p, draw_edge(N[i], N[j], string("x<sub>",i,j,"</sub>")))
end
end
compose(context(), fill(nothing), stroke("black"), p)
end
Out[2]:
One standard pooling test instance is the Haverly model,
In [3]:
N = PoolingNetwork(3, 1, 2,
sparse([1,2,3,3,4,4],[4,4,5,6,5,6],[6,16,1,5,-9,-15],6,6),
[3,1,2,0,2.5,1.5],
[1,1,1,1,1/3,2/3] );
with quality parameters
In [4]:
N.q
Out[4]:
capacities,
In [5]:
N.cap
Out[5]:
and costs
In [6]:
N.cost
Out[6]:
Edges to terminals have a negative cost, corresponding to an income when a commodity is delivered.
In [7]:
draw_network(N)
Out[7]:
In [8]:
using Polyopt
In [9]:
x14, x24, x35, x36, x45, x46, w4 = variables(["x14"; "x24"; "x35"; "x36"; "x45"; "x46"; "w4"]);
In the minimium cost formulation we minimize $\sum_{(i,j)\in N} c_{i,j}x_{i,j}$.
In [10]:
obj = N.cost[1,4]*x14 + N.cost[2,4]*x24 + N.cost[3,5]*x35 + N.cost[3,6]*x36 + N.cost[4,5]*x45 + N.cost[4,6]*x46;
We then have flow equilibrium constraints $x_{14} + x_{24} = x_{45} + x_{46}$,
In [11]:
flow_eq = x14 + x24 - (x45 + x46);
and a (nonconvex) defining equation for the blending variables $w_4=(x_{45} + x_{46}) = q_1x_{14} + q_2 x_{24}$ as the blended quality at the pool,
In [12]:
blend_eq = w4*(x45 + x46) - (N.q[1]*x14 + N.q[2]*x24);
quality bounds at the terminals $q_5(x_{45}+x_{36}) \geq (w_4 x_{45} + q_3 x_{35})$ and $q_6(x_{46}+x_{36}) \geq (w_4 x_{46} + q_3 x_{36})$
In [13]:
qual_bnd = [ N.q[5]*(x45+x36)-(w4*x45+N.q[3]*x35); N.q[6]*(x46+x36)-(w4*x46+N.q[3]*x36) ];
capacity bounds,
In [14]:
cap_bnd = [ N.cap[1]-x14; N.cap[2]-x24; N.cap[3]-(x35+x36); N.cap[4]-(x45+x46); N.cap[5]-(x35+x45); N.cap[6]-(x36+x46) ];
and finally nonnegativity of flow variables,
In [15]:
flow_bnd = [ x14; x24; x35; x36; x45; x46 ];
We solve the second-order moment relaxation
In [16]:
prob = momentprob(2, obj, [qual_bnd; cap_bnd; flow_bnd], [flow_eq; blend_eq]);
In [17]:
X, Z, t, y, solsta = solve_mosek(prob);
We examine the solution,
In [18]:
[prob.basis y]
Out[18]:
and we we test feasibility,
In [19]:
xo = Polyopt.vectorize([x14,x24,x35,x36,x45,x46,w4],2*2)*y
Out[19]:
In [20]:
[ Polyopt.evalpoly(h, xo) for h=[flow_eq; blend_eq] ]
Out[20]:
In [21]:
[ Polyopt.evalpoly(g, xo) for g=[qual_bnd; cap_bnd; flow_bnd] ]
Out[21]:
since the extracted monomial solution is feasible, we have found the optimal solution.