Add missing problems

data/tetra.jl

@@ -0,0 +1,25 @@
+xe_tetra = [
+    0 0 0
+    1 0 0
+    0 1 0
+    0 0 1
+    5//10 2//10 1//10
+];
+
+Tets_tetra = [
+    1 2 3 5
+    5 4 3 1
+    5 4 1 2
+    5 2 3 4
+];
+
+Tets_tetra = vec(reshape(Tets_tetra, 16, 1));
+xe_tetra = vec(reshape(xe_tetra, 15, 1));
+Constants_tetra = [1, 2, 3, 4];
+
+include("tetra_duct12.jl")
+include("tetra_duct15.jl")
+include("tetra_duct20.jl")
+include("tetra_hook.jl")
+include("tetra_foam5.jl")
+include("tetra_gear.jl")

data/triangle.jl

@@ -0,0 +1,20 @@
+xe = [
+    0 0
+    1 0
+    0 1
+    64//100 25//100
+];
+
+Tr = [
+    1 2 4
+    4 2 3
+    4 3 1
+];
+
+Tr = vec(reshape(Tr, 9, 1));
+xe = vec(reshape(xe, 8, 1));
+Constants = [1, 2, 3];
+
+include("triangle_deer.jl")
+include("triangle_pacman.jl")
+include("triangle_turtle.jl")

src/ExaModelsExamples.jl

@@ -15,13 +15,14 @@ include("quadrotor.jl")
 include("goddard.jl")
 include("robot.jl")
 include("rocket.jl")
-include("bearing.jl") 
+include("bearing.jl")
 include("camshape.jl")
 include("elec.jl")
 include("steering.jl")
 include("pinene.jl")
 include("marine.jl")
 include("gasoil.jl")
+include("pde_models.jl")
 
 const NAMES = filter(names(ExaModelsExamples; all = true)) do x
     str = string(x)

src/catmix.jl

@@ -0,0 +1,74 @@
+# Catalyst Mixing Problem
+# Collocation formulation
+# COPS 3.0 - November 2002
+# COPS 3.1 - March 2004
+
+function catmix_model(nh; 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
+
+

src/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 chain_model(n; 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
+

src/glider.jl

@@ -0,0 +1,100 @@
+# 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; T = Float64, backend = nothing, kwargs...)
+    # 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
+
+    c = ExaModels.ExaCore(T; backend = backend)
+    t_f = ExaModels.variable(c, 1; lvar = 0, start = 1)
+    x = ExaModels.variable(c, nh+1; lvar = zeros(nh+1), start = [x_0 + vx_0*(k/nh) for k in 0:nh])
+    y = ExaModels.variable(c, nh+1; start = [y_0 + (k/nh)*(y_f - y_0) for k in 0:nh])
+    vx = ExaModels.variable(c, nh+1; lvar = zeros(nh+1), start = fill(vx_0, nh+1))
+    vy = ExaModels.variable(c, nh+1; start = fill(vy_0, nh+1))
+    cL = ExaModels.variable(c, nh+1; lvar = fill(cL_min,nh+1), uvar = fill(cL_max, nh+1), start = fill(cL_max/2, nh+1))
+
+    ExaModels.objective(c, -x[nh+1])
+
+    ExaModels.constraint(
+        c,
+        x[j] - (x[j-1] + 0.5/nh * (vx[j] + vx[j-1]) * t_f[1]) for j in 2:nh+1
+    )
+
+    ExaModels.constraint(
+        c,
+        y[j] - (y[j-1] + 0.5 * t_f[1]/nh * (vy[j] + vy[j-1])) for j in 2:nh+1
+    )
+
+    ExaModels.constraint(
+        c,
+        vx[j] - (vx[j-1] + 0.5 * t_f[1]/nh * (vx_dot[j] + vx_dot[j-1])) for j in 2:nh+1
+    )
+
+    ExaModels.constraint(
+        c,
+        vy[j] - (vy[j-1] + 0.5 * t_f[1]/nh * (vy_dot[j] + vy_dot[j-1])) for j in 2:nh+1
+    )
+
+    ExaModels.constraint(
+        c,
+        x[1] - x_0
+    )
+
+    ExaModels.constraint(
+        c,
+        y[1] - y_0
+    )
+
+    ExaModels.constraint(
+        c,
+        y[nh+1] - y_f
+    )
+
+    ExaModels.constraint(
+        c,
+        vx[1] - vx_0
+    )
+
+    ExaModels.constraint(
+        c,
+        vx[nh+1] - vx_f
+    )
+
+    ExaModels.constraint(
+        c,
+        vy[1] - vy_0
+    )
+
+    ExaModels.constraint(
+        c,
+        vy[nh+1] - vy_f
+    )
+
+    ExaModels.Examodel(c,; kwargs...)
+end
+
+
+println("Running tests...")
+using NLPModelsIpopt, ExaModels, MadNLP
+
+result = madnlp((glider_model(200));print_level = MadNLP.INFO )