From 673c027d3a4508caf81f778fd809111f80769609 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Pacaud?= Date: Thu, 22 Feb 2024 12:04:14 +0100 Subject: [PATCH] add channel, glider and methanol instances (#3) --- src/COPSBenchmark.jl | 3 + src/channel.jl | 81 +++++++++++++++++++++++++ src/glider.jl | 72 ++++++++++++++++++++++ src/methanol.jl | 139 +++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 3 + 5 files changed, 298 insertions(+) create mode 100644 src/channel.jl create mode 100644 src/glider.jl create mode 100644 src/methanol.jl diff --git a/src/COPSBenchmark.jl b/src/COPSBenchmark.jl index 61253c6..479341e 100644 --- a/src/COPSBenchmark.jl +++ b/src/COPSBenchmark.jl @@ -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") diff --git a/src/channel.jl b/src/channel.jl new file mode 100644 index 0000000..0ce8a53 --- /dev/null +++ b/src/channel.jl @@ -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 + diff --git a/src/glider.jl b/src/glider.jl new file mode 100644 index 0000000..6fe6f4b --- /dev/null +++ b/src/glider.jl @@ -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 + diff --git a/src/methanol.jl b/src/methanol.jl new file mode 100644 index 0000000..07b1e44 --- /dev/null +++ b/src/methanol.jl @@ -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 + diff --git a/test/runtests.jl b/test/runtests.jl index b78931a..b41a12a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -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),