add COPS instances formulated with ExaModels (#11)

Author: frapac

Project.toml

@@ -1,12 +1,20 @@
 name = "COPSBenchmark"
 uuid = "3c263af4-7bba-497b-afe4-00c59a0a48a5"
-authors = ["fpacaud <[email protected]>"]
-version = "0.1.0"
+version = "0.2.0"
+authors = ["MadSuite <https://madsuite.org/>"]
 
 [deps]
-JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
+[weakdeps]
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2"
+
+[extensions]
+COPSBenchmarkExaModels = "ExaModels"
+COPSBenchmarkJuMP = "JuMP"
+
 [compat]
+ExaModels = "0.9"
 JuMP = "^1.19"
-julia = "1.6"
+julia = "1.10"

ext/COPSBenchmarkExaModels/COPSBenchmarkExaModels.jl

@@ -0,0 +1,25 @@
+module COPSBenchmarkExaModels
+
+using Random
+import COPSBenchmark
+import COPSBenchmark: ExaModelsBackend
+using ExaModels
+
+include("bearing.jl")
+include("camshape.jl")
+include("catmix.jl")
+include("chain.jl")
+include("elec.jl")
+include("gasoil.jl")
+include("glider.jl")
+include("marine.jl")
+include("methanol.jl")
+include("minsurf.jl")
+include("pinene.jl")
+include("robot.jl")
+include("rocket.jl")
+include("steering.jl")
+include("pde_models.jl")
+
+end
+

ext/COPSBenchmarkExaModels/bearing.jl

@@ -0,0 +1,44 @@
+# Journal bearing problem
+# Michael Merritt - Summer 2000
+# COPS 2.0 - September 2000
+# COPS 3.0 - November 2002
+# COPS 3.1 - March 2004
+
+function COPSBenchmark.bearing_model(nx, ny, ::ExaModelsBackend; T = Float64, backend = nothing, kwargs...)
+    b = 10              # grid is (0,2*pi)x(0,2*b)
+    e = 0.1             # eccentricity
+
+    hx = 2*pi / (nx+1)  # grid spacing
+    hy = 2*b / (ny+1)   # grid spacing
+    area = 0.5*hx*hy    # area of triangle
+
+    wq(i) = (1.0 + e*cos((i-1)*hx))^3
+    v0 = [max(sin((i-1)*hx), 0.0) for i in 1:nx+2, j in 1:ny+2]
+
+    core = ExaModels.ExaCore(T; backend=backend)
+
+    v = ExaModels.variable(core, 1:nx+2, 1:ny+2; lvar = 0.0, start=v0)
+
+    ExaModels.objective(
+        core,
+        0.5*(hx*hy/6.0) * (wq(i) + 2*wq(i+1))*(((v[i+1,j]-v[i,j])/hx)^2 + ((v[i,j+1]-v[i,j])/hy)^2) for i in 1:nx+1, j in 1:ny+1
+    )
+    ExaModels.objective(
+        core,
+        0.5*(hx*hy/6.0) * (2*wq(i) + 2*wq(i-1))*(((v[i-1,j]-v[i,j])/hx)^2 + ((v[i,j-1]-v[i,j])/hy)^2) for i in 2:nx+2, j in 2:ny+2
+    )
+    ExaModels.objective(
+        core,
+        -hx*hy*e*sin((i-1)*hx)*v[i, j] for i in 1:nx+2, j in 1:ny+2
+    )
+
+    ExaModels.constraint(core, v[i, 1] for i in 1:nx+2)
+    ExaModels.constraint(core, v[i, ny+2] for i in 1:nx+2)
+    ExaModels.constraint(core, v[1, i] for i in 1:ny+2)
+    ExaModels.constraint(core, v[nx+2, i] for i in 1:ny+2)
+
+    return ExaModels.ExaModel(core; kwargs...)
+end
+
+
+

ext/COPSBenchmarkExaModels/camshape.jl

@@ -0,0 +1,60 @@
+# Cam Shape Problem
+# Alexander S. Bondarenko - Summer 1998
+# COPS 2.0 - September 2000
+# COPS 3.0 - November 2002
+# COPS 3.1 - March 2004
+
+function COPSBenchmark.camshape_model(n, ::ExaModelsBackend; T = Float64, backend = nothing, kwargs...)
+
+    R_v = 1.0         # design parameter related to the valve shape
+    R_max = 2.0       # maximum allowed radius of the cam
+    R_min = 1.0       # minimum allowed radius of the cam
+    alpha = 1.5       # curvature limit parameter
+
+    d_theta = 2*pi/(5*(n+1))   # angle between discretization points
+
+    core = ExaModels.ExaCore(T; backend= backend, minimize=false)
+
+    # radius of the cam at discretization points
+    r= ExaModels.variable(core, 1:n; lvar =R_min, uvar = R_max, start=(R_min+R_max)/2.0)
+
+    ExaModels.objective(core, (pi*R_v)/n * r[i] for i in 1:n)
+
+    # Convexity
+    ExaModels.constraint(
+        core,
+        - r[i-1]*r[i] - r[i]*r[i+1] + 2*r[i-1]*r[i+1]*cos(d_theta) for i=2:n-1; lcon = -Inf, ucon = 0.0,
+            )
+    ExaModels.constraint(
+        core,
+        - R_min*r[1] - r[1]*r[2] + 2*R_min*r[2]*cos(d_theta); lcon = -Inf, ucon = 0.0
+    )
+    ExaModels.constraint(
+        core,
+        - R_min^2 - R_min*r[1] + 2*R_min*r[1]*cos(d_theta); lcon = -Inf, ucon = 0.0
+    )
+    ExaModels.constraint(
+        core,
+        - r[n-1]*r[n] - r[n]*R_max + 2*r[n-1]*R_max*cos(d_theta); lcon = -Inf, ucon = 0.0
+    )
+    ExaModels.constraint(
+        core,
+        - 2*R_max*r[n] + 2*r[n]^2*cos(d_theta); lcon = -Inf, ucon = 0.0
+    )
+    # Curvature
+    ExaModels.constraint(
+        core,
+        (r[i+1] - r[i]) for i=1:n-1; lcon = -alpha*d_theta, ucon = alpha*d_theta,
+    )
+    ExaModels.constraint(
+        core,
+        (r[1] - R_min); lcon = -alpha*d_theta, ucon = alpha*d_theta
+    )
+    ExaModels.constraint(
+        core,
+        (R_max - r[n]); lcon = -alpha*d_theta, ucon = alpha*d_theta
+    )
+
+    return ExaModels.ExaModel(core; kwargs...)
+end
+

ext/COPSBenchmarkExaModels/catmix.jl

@@ -0,0 +1,74 @@
+# Catalyst Mixing Problem
+# Collocation formulation
+# COPS 3.0 - November 2002
+# COPS 3.1 - March 2004
+
+function COPSBenchmark.catmix_model(nh, ::ExaModelsBackend; T = Float64, backend = nothing, kwargs...)
+    ne = 2
+    nc = 3
+
+    tf = 1
+    h = tf / nh   # Final time
+
+    rho = [
+        0.11270166537926,
+        0.50000000000000,
+        0.88729833462074,
+    ]
+    bc = [1.0, 0.0]  # Boundary conditions for x
+    alpha = 0.0      # Smoothing parameter
+    rho_index = [(i, rho[i]) for i in 1:nc]
+
+    c = ExaModels.ExaCore(T; backend = backend)
+    u = ExaModels.variable(c, nh, nc; lvar = zeros(nh, nc), uvar = ones(nh, nc), start = zeros(nh, nc))
+    v = ExaModels.variable(c, nh, ne; start = [mod(j, ne) for i in 1:nh, j in 1:ne])
+    w = ExaModels.variable(c, nh, nc, ne; start = zeros(nh, nc, ne))
+    pp = ExaModels.variable(c, nh, nc, ne; start = [mod(k, ne) for i in 1:nh, j in 1:nc, k in 1:ne])
+    Dpp = ExaModels.variable(c, nh, nc, ne; start = zeros(nh, nc, ne))
+    ppf = ExaModels.variable(c, ne; start = [mod(i,ne) for i in 1:ne])
+
+    ExaModels.objective(c, -1.0 + ppf[1] + ppf[2])
+    ExaModels.objective(c, alpha/h*(u[i+1, j] - u[i, j])^2 for i in 1:nh-1, j in 1:nc)
+
+    ExaModels.constraint(
+        c,
+        pp[i, k, s] - v[i, s] - h*sum(w[i, j, s]*(rho^j/factorial(j)) for j in 1:nc) for i=1:nh,  (k, rho) in rho_index, s=1:ne
+    )
+
+    ExaModels.constraint(
+        c,
+        Dpp[i, k, s] - sum(w[i, j, s]*(rho^(j-1)/factorial(j-1)) for j in 1:nc) for i=1:nh, (k, rho) in rho_index, s=1:ne
+    )
+
+    ExaModels.constraint(
+        c,
+        ppf[s] - v[nh, s] - h * sum(w[nh, j, s] / factorial(j) for j in 1:nc) for s in 1:ne
+    )
+
+    ExaModels.constraint(
+        c,
+        v[i, s] + sum(w[i, j, s] * h / factorial(j) for j in 1:nc) - v[i+1, s] for i in 1:nh-1, s in 1:ne
+    )
+
+
+
+    ExaModels.constraint(
+        c,
+        Dpp[i,j,1] - u[i,j] * (10.0*pp[i,j,2] - pp[i,j,1]) for i=1:nh, j=1:nc
+    )
+
+    ExaModels.constraint(
+        c,
+        Dpp[i,j,2] - u[i,j] * (pp[i,j,1] - 10.0*pp[i,j,2]) + (1 - u[i,j])*pp[i,j,2] for i=1:nh, j=1:nc
+    )
+
+
+    ExaModels.constraint(
+        c,
+        v[1, s] - bc for (s, bc) in [(i, bc[i]) for i in 1:ne]
+    )
+
+    return ExaModels.ExaModel(c; kwargs...)
+end
+
+

ext/COPSBenchmarkExaModels/chain.jl

@@ -0,0 +1,72 @@
+# Hanging Chain
+
+# Find the chain (of uniform density) of length L suspended between two points with minimal
+# potential energy.
+
+#   This is problem 4 in the COPS (Version 3) collection of
+#   E. Dolan and J. More'
+#   see "Benchmarking Optimization Software with COPS"
+#   Argonne National Labs Technical Report ANL/MCS-246 (2004)
+
+function COPSBenchmark.chain_model(n, ::ExaModelsBackend; T = Float64, backend = nothing, kwargs...)
+    nh = max(2, div(n - 4, 4))
+
+    L = 4
+    a = 1
+    b = 3
+    tmin = b > a ? 1 / 4 : 3 / 4
+    tf = 1.0
+    h = tf / nh
+
+    c = ExaModels.ExaCore(T; backend = backend)
+    u = ExaModels.variable(c, nh + 1; start = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1])
+    x1 = ExaModels.variable(c, nh + 1; start = [4 * abs(b - a) * k / nh * (1 / 2 * k / nh - tmin) + a for k in 1:nh+1])
+    x2 = ExaModels.variable(c, nh + 1; start = [(4 * abs(b - a) * k / nh * (1 / 2 * k / nh - tmin) + a) *
+        (4 * abs(b - a) * (k / nh - tmin)) for k in 1:nh+1])
+    x3 = ExaModels.variable(c, nh + 1;  start = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1])
+
+    ExaModels.objective(c, x2[nh + 1])
+
+    ExaModels.constraint(
+        c,
+        x1[j + 1] - x1[j] - 1 / 2 * h * (u[j] + u[j + 1]) for j in 1:nh
+    )
+
+    ExaModels.constraint(
+        c,
+        x1[1] - a
+    )
+
+    ExaModels.constraint(
+        c,
+        x1[nh + 1] - b
+    )
+
+    ExaModels.constraint(
+        c,
+        x2[1]
+    )
+
+    ExaModels.constraint(
+        c,
+        x3[1]
+    )
+
+    ExaModels.constraint(
+        c,
+        x3[nh+1] - L
+    )
+
+    ExaModels.constraint(
+        c,
+        x2[j + 1] - x2[j] - 1 / 2 * h * (x1[j] * sqrt(1 + u[j]^2) + x1[j + 1] * sqrt(1 + u[j + 1]^2)) for j in 1:nh
+    )
+
+    ExaModels.constraint(
+        c,
+        x3[j + 1] - x3[j] - 1 / 2 * h * (sqrt(1 + u[j]^2) + sqrt(1 + u[j + 1]^2)) for j in 1:nh
+    )
+
+    return ExaModels.ExaModel(c; kwargs...)
+end
+

ext/COPSBenchmarkExaModels/channel.jl.bkp

@@ -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; T = Float64, backend = nothing, kwargs...)
+    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
+
+    core = ExaModels.ExaCore(T; backend= backend)
+
+    v = ExaModels.variable(model, 1:nh, 1:nd)
+    w = ExaModels.variable(model, 1:nh, 1:nc; start=0.0)
+
+    uc = ExaModels.variable(model, 1:nh, 1:nc, 1:nd; start=[v0[i, s] for i=1:nh, j=1:nc, s=1:nd])
+    Duv = ExaModels.variable(model, 1:nh, 1:nc, 1:nd; start=0.0)
+    
+    # Constant objective
+    ExaModels.objective(model, Min, 1.0)
+
+    # Collocation model
+    ExaModels.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),
+    )
+    ExaModels.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
+    ExaModels.constraint(model, bc_1, v[1, 1] == bc[1, 1])
+    ExaModels.constraint(model, bc_2, v[1, 2] == bc[2, 1])
+    ExaModels.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]
+    )
+    ExaModels.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]
+    )
+    ExaModels.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]
+    )
+    ExaModels.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 ExaModels.ExaModel(core; kwargs...)
+end
+