Skip to content

Commit

Permalink
add channel, glider and methanol instances (#3)
Browse files Browse the repository at this point in the history
  • Loading branch information
frapac authored Feb 22, 2024
1 parent 15df78b commit 673c027
Show file tree
Hide file tree
Showing 5 changed files with 298 additions and 0 deletions.
3 changes: 3 additions & 0 deletions src/COPSBenchmark.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,12 @@ using JuMP
include("bearing.jl")
include("camshape.jl")
include("catmix.jl")
include("channel.jl")
include("elec.jl")
include("gasoil.jl")
include("glider.jl")
include("marine.jl")
include("methanol.jl")
include("pinene.jl")
include("robot.jl")
include("rocket.jl")
Expand Down
81 changes: 81 additions & 0 deletions src/channel.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,81 @@
# Flow in a Channel Problem
# Collocation formulation
# Alexander S. Bondarenko - Summer 1998
# COPS 2.0 - September 2000
# COPS 3.0 - November 2002
# COPS 3.1 - March 2004

function channel_model(nh)
nc = 4
nd = 4
R = 10.0 # Reynolds number
tf = 1.0
h = tf / nh

bc = [0.0 1.0; 0.0 0.0]
rho = [0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703]
t = [(i-1)*h for i in 1:nh+1]

# Initial value
v0 = zeros(nh, nd)
for i in 1:nh
v0[i, 1] = t[i]^2*(3.0 - 2.0*t[i])
v0[i, 2] = 6*t[i]*(1.0 - t[i])
v0[i, 3] = 6*(1.0 - 2.0*t[i])
v0[i, 4] = -12.0
end

model = Model()

@variable(model, v[i=1:nh, j=1:nd])
@variable(model, w[i=1:nh, j=1:nc], start=0.0)

@variable(model, uc[i=1:nh, j=1:nc, s=1:nd], start=v0[i, s])
@variable(model, Duc[i=1:nh, j=1:nc, s=1:nd], start=0.0)

# Constant objective
@objective(model, Min, 1.0)

# Collocation model
@constraint(
model,
[i=1:nh, j=1:nc, s=1:nd],
uc[i, j, s] == v[i,s] + h*sum(w[i,k]*(rho[j]^k/factorial(k)) for k in 1:nc),
)
@constraint(
model,
[i=1:nh, j=1:nc, s=1:nd],
Duc[i, j, s] == sum(v[i,k]*((rho[j]*h)^(k-s)/factorial(k-s)) for k in s:nd) +
h^(nd-s+1) * sum(w[i, k]*(rho[j]^(k+nd-s)/factorial(k+nd-s)) for k in 1:nc)
)
# Boundary
@constraint(model, bc_1, v[1, 1] == bc[1, 1])
@constraint(model, bc_2, v[1, 2] == bc[2, 1])
@constraint(
model,
bc_3,
sum(v[nh, k]*(h^(k-1)/factorial(k-1)) for k in 1:nd) +
h^nd * sum(w[nh, k]/factorial(k+nd-1) for k in 1:nc) == bc[1, 2]
)
@constraint(
model,
bc_4,
sum(v[nh, k]*(h^(k-2)/factorial(k-2)) for k in 2:nd) +
h^(nd-1) * sum(w[nh, k]/factorial(k+nd-2) for k in 1:nc) == bc[2, 2]
)
@constraint(
model,
continuity[i=1:nh-1, s=1:nd],
sum(v[i, k]*(h^(k-s)/factorial(k-s)) for k in s:nd)
+ h^(nd-s+1)* sum(w[i, k]/factorial(k+nd-s) for k in 1:nc) == v[i+1, s]
)
@constraint(
model,
collocation[i=1:nh, j=1:nc],
sum(w[i, k] * (rho[j]^(k-1)/factorial(k-1)) for k in 1:nc) ==
R * (Duc[i, j, 2] * Duc[i, j, 3] - Duc[i, j, 1] * Duc[i, j, 4])
)

return model
end

72 changes: 72 additions & 0 deletions src/glider.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@
# Hang Glider Problem
# Trapezoidal formulation
# David Bortz - Summer 1998
# COPS 2.0 - September 2000
# COPS 3.0 - November 2002
# COPS 3.1 - March 2004

function glider_model(nh)
# Design parameters
x_0 = 0.0
y_0 = 1000.0
y_f = 900.0
vx_0 = 13.23
vx_f = 13.23
vy_0 = -1.288
vy_f = -1.288
u_c = 2.5
r_0 = 100.0
m = 100.0
g = 9.81
c0 = 0.034
c1 = 0.069662
S = 14.0
rho = 1.13
cL_min = 0.0
cL_max = 1.4

model = Model()

@variable(model, step >= 0.0, start=1/nh)
# State variables
@variable(model, 0.0 <= x[k=0:nh], start=x_0 + vx_0*(k/nh))
@variable(model, y[k=0:nh], start=y_0 + (k/nh)*(y_f - y_0))
@variable(model, 0.0 <= vx[k=0:nh], start=vx_0)
@variable(model, vy[k=0:nh], start=vy_0)
# Control variables
@variable(model, cL_min <= cL[k=0:nh] <= cL_max, start=cL_max/2.0)

@objective(model, Max, x[nh])

@expressions(model, begin
r[i=0:nh], (x[i]/r_0 - 2.5)^2
u[i=0:nh], u_c*(1 - r[i])*exp(-r[i])
w[i=0:nh], vy[i] - u[i]
v[i=0:nh], sqrt(vx[i]^2 + w[i]^2)
D[i=0:nh], 0.5*(c0+c1*cL[i]^2)*rho*S*v[i]^2
L[i=0:nh], 0.5*cL[i]*rho*S*v[i]^2
vx_dot[i=0:nh], (-L[i]*(w[i]/v[i]) - D[i]*(vx[i]/v[i]))/m
vy_dot[i=0:nh], (L[i]*(vx[i]/v[i]) - D[i]*(w[i]/v[i]))/m - g
end)

# Dynamics
@constraints(model, begin
x_eqn[j=1:nh], x[j] == x[j-1] + 0.5 * step * (vx[j] + vx[j-1])
y_eqn[j=1:nh], y[j] == y[j-1] + 0.5 * step * (vy[j] + vy[j-1])
vx_eqn[j=1:nh], vx[j] == vx[j-1] + 0.5 * step * (vx_dot[j] + vx_dot[j-1])
vy_eqn[j=1:nh], vy[j] == vy[j-1] + 0.5 * step * (vy_dot[j] + vy_dot[j-1])
end)
# Boundary constraints
@constraints(model, begin
x[0] == x_0
y[0] == y_0
y[nh] == y_f
vx[0] == vx_0
vx[nh] == vx_f
vy[0] == vy_0
vy[nh] == vy_f
end)

return model
end

139 changes: 139 additions & 0 deletions src/methanol.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,139 @@
# Methanol-to-Hydrocarbons Problem
# Collocation formulation
# Michael Merritt - Summer 2000
# COPS 2.0 - September 2000
# COPS 3.0 - November 2002
# COPS 3.1 - March 2004

function methanol_model(nh)
ne = 3
np = 5
nc = 3
nm = 17

rho = [0.11270166537926, 0.5, 0.88729833462074]
# times at which observations made
tau = [
0,
0.050,
0.065,
0.080,
0.123,
0.233,
0.273,
0.354,
0.397,
0.418,
0.502,
0.553,
0.681,
0.750,
0.916,
0.937,
1.122,
]
tf = tau[nm] # ODEs defined in [0,tf]
h = tf / nh # uniform interval length
t = [(i-1)*h for i in 1:nh+1] # partition

# itau[i] is the largest integer k with t[k] <= tau[i]
itau = Int[min(nh, floor(tau[i]/h)+1) for i in 1:nm]

# Concentrations
z = [
1.0000 0 0;
0.7085 0.1621 0.0811;
0.5971 0.1855 0.0965;
0.5537 0.1989 0.1198;
0.3684 0.2845 0.1535;
0.1712 0.3491 0.2097;
0.1198 0.3098 0.2628;
0.0747 0.3576 0.2467;
0.0529 0.3347 0.2884;
0.0415 0.3388 0.2757;
0.0261 0.3557 0.3167;
0.0208 0.3483 0.2954;
0.0085 0.3836 0.2950;
0.0053 0.3611 0.2937;
0.0019 0.3609 0.2831;
0.0018 0.3485 0.2846;
0.0006 0.3698 0.2899;
]

bc = [1.0, 0.0, 0.0]

# Starting-value
v0 = zeros(nh, ne)
for i in 1:itau[1], s in 1:ne
v0[i, s] = bc[s]
end
for j in 2:nm, i in itau[j-1]+1:itau[j], s in 1:ne
v0[i, s] = z[j, s]
end
for i in itau[nm]+1:nh, s in 1:ne
v0[i, s] = z[nm, s]
end
v0 .= 0.001

model = Model()

@variable(model, theta[1:np] >= 0.0, start=1.0) # ODE parameters
@variable(model, v[i=1:nh, s=1:ne], start=v0[i, s])
@variable(model, w[i=1:nh, j=1:nc, s=1:ne], start=0.0)
@variable(model, uc[i=1:nh, j=1:nc, s=1:ne], start=v0[i, s])
@variable(model, Duc[i=1:nh, j=1:nc, s=1:ne], start=0.0)

# error
@expression(
model,
error[j=1:nm, s=1:ne],
v[itau[j],s] + sum(w[itau[j],k,s]*(tau[j]-t[itau[j]])^k/(factorial(k)*h^(k-1)) for k in 1:nc) - z[j,s]
)

# L2 error
@objective(model, Min, sum(error[j, s]^2 for j in 1:nm, s in 1:ne))

# Collocation model
@constraint(
model,
[i=1:nh, j=1:nc, s=1:ne],
uc[i, j, s] == v[i,s] + h*sum(w[i,k,s]*(rho[j]^k/factorial(k)) for k in 1:nc),
)
@constraint(
model,
[i=1:nh, j=1:nc, s=1:ne],
Duc[i, j, s] == sum(w[i,k,s]*(rho[j]^(k-1)/factorial(k-1)) for k in 1:nc),
)
# Boundary
@constraint(model, bc_cond[s=1:ne], v[1, s] == bc[s])
# Continuity
@constraint(
model,
[i=1:nh-1, s=1:ne],
v[i, s] + sum(w[i, j, s]*h/factorial(j) for j in 1:nc) == v[i+1, s]
)
# Dynamics
@constraint(
model,
collocation_eqn1[i=1:nh, j=1:nc],
Duc[i,j,1] == - (2*theta[2] - (theta[1]*uc[i,j,2])/((theta[2]+theta[5])*uc[i,j,1]+uc[i,j,2]) +
theta[3] + theta[4])*uc[i,j,1]
)
@constraint(
model,
collocation_eqn2[i=1:nh, j=1:nc],
Duc[i,j,2] == (theta[1]*uc[i,j,1]*(theta[2]*uc[i,j,1]-uc[i,j,2]))/
((theta[2]+theta[5])*uc[i,j,1]+uc[i,j,2]) +
theta[3]*uc[i,j,1]
)
@constraint(
model,
collocation_eqn3[i=1:nh, j=1:nc],
Duc[i,j,3] == (theta[1]*uc[i,j,1]*(uc[i,j,2]+theta[5]*uc[i,j,1]))/
((theta[2]+theta[5])*uc[i,j,1]+uc[i,j,2]) +
theta[4]*uc[i,j,1]
)

return model
end

3 changes: 3 additions & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,9 +8,12 @@ COPS_INSTANCES = [
(COPSBenchmark.bearing_model, (50, 50), -1.5482e-1),
(COPSBenchmark.camshape_model, (1000,), 4.2791), # TODO: result is slightly different
(COPSBenchmark.catmix_model, (100,), -4.80556e-2),
(COPSBenchmark.channel_model, (200,), 1.0),
(COPSBenchmark.elec_model, (50,), 1.0552e3),
(COPSBenchmark.gasoil_model, (100,), 5.2366e-3),
(COPSBenchmark.glider_model, (100,), 1.25505e3),
(COPSBenchmark.marine_model, (100,), 1.97462e7),
(COPSBenchmark.methanol_model, (100,), 9.02229e-3),
(COPSBenchmark.pinene_model, (100,), 1.98721e1),
(COPSBenchmark.robot_model, (200,), 9.14138),
(COPSBenchmark.rocket_model, (400,), 1.01283),
Expand Down

0 comments on commit 673c027

Please sign in to comment.