Add-goddardrocket

src/ADNLPProblems/minimalsurface.jl

@@ -0,0 +1,105 @@
+using Plots
+using ADNLPModels, NLPModels, NLPModelsIpopt, DataFrames, LinearAlgebra, Distances, SolverCore, PyPlot
+
+function minimalsurface(; n::Int = default_nvar, type::Val{T} = Val(Float64), kwargs...) where {T}
+
+    # domain definition
+    xmin = T(0.)
+    xmax = T(1.)
+    ymin = T(0.)
+    ymax = T(1.)
+
+    # Definition of the mesh
+    nx = 20 # number of points according to the direction x
+    ny = 20 # number of points according to the direction y
+
+
+    x_mesh = LinRange(xmin, xmax, nx) # coordinates of the mesh points x
+    y_mesh = LinRange(ymin, ymax, ny) # coordinates of the mesh points y
+
+    v_D = zeros(nx,ny) # Surface matrix initialization
+    for i in 1:nx
+        for j in 1:ny
+            v_D[i, j] = T(1 - (2 * x_mesh[i] - 1)^2)
+        end 
+    end  
+
+
+    function Objective(v)
+        v_reshape = reshape(v, (nx, ny)) # vector to matrix conversion
+        hx = T(1/nx)  # step on the x axis
+        hy = T(1/ny)  # step on the y axis
+        S = T(0.) # sum initialization
+        # Calculation of the gradient according to x
+        for i in 1:nx
+            for j in 1:ny
+                if i == 1
+                    gi = T((v_reshape[i+1, j] - v_reshape[i, j])/hx)
+                elseif i == nx
+                    gi = T((v_reshape[i, j] - v_reshape[i-1, j])/hx)
+                else 
+                    gi = T((v_reshape[i+1, j] - v_reshape[i, j])/(2 * hx))
+                end 
+        # Calculation of the gradient according to x
+                if j == 1
+                    gj = T((v_reshape[i, j+1] - v_reshape[i, j])/hy)            
+                elseif j == ny
+                    gj = T((v_reshape[i, j] - v_reshape[i, j-1])/hy) 
+                else 
+                    gj = T((v_reshape[i, j+1] - v_reshape[i, j])/(2 * hy)) 
+                end 
+
+                s = T(hx * hy * (sqrt(1 + (gi^2) +(gj)^2))) # Approximation of the derivative with the method of rectangles
+                S = S + s # Update the value of S
+            end
+        end
+        return(S)
+    end 
+
+    function constraints(v)
+        v_reshape = reshape(v, (nx, ny)) # vector to matrix conversion
+        c = similar(v_reshape, nx*ny + 2*(nx +ny)) # creating a constraint vector
+        index = 1
+        v_L = zeros(T, nx,ny) # Creation of an obstacle called v_L
+        for i in 1:nx
+            for j in 1:ny
+                if norm(x_mesh[i]-(1/2)) ≤ 1/4 && norm(y_mesh[j]-(1/2)) ≤ 1/4
+                    v_L[i, j] = T(1.) # Update the value of v_L according to the values ​​of x and y
+                end
+            end 
+        end  
+        for i in 1:nx
+            for j in 1:ny
+                c[index] = T(v_reshape[i, j] - v_L[i, j]) # Constraint that the surface must be above the obstruction
+                index = index + 1
+            end 
+        end 
+        for j in 1:ny
+            c[index] = T(v_reshape[1, j]) # Constraint: when x=1 or x=nx, the surface is set to 0
+            index = index + 1
+            c[index] =  T(v_reshape[nx, j]) # Constraint: when x=1 or x=nx, the surface is set to 0
+            index = index + 1
+        end 
+        for i in 1:nx
+            c[index] = T(v_reshape[i, 1] - 1 + (2 * i -1)^2) # Constraint: when y=1 or y=ny, the surface follows the function " 1 + (2 * x -1)^2 "
+            index = index + 1
+            c[index] = T(v_reshape[i, ny] - 1 + (2 * i -1)^2) # Constraint: when y=1 or y=ny, the surface follows the function " 1 + (2 * x -1)^2 "
+            index = index + 1
+
+        end 
+        return c
+    end
+
+
+    lcon = zeros(T, nx * ny + 2 * nx + 2 * ny) # Lower bound all equal to 0
+    ucon = zeros(T, nx * ny + 2 * nx + 2 * ny) # Creation of the upper bound vector
+    ucon[1 : ny * nx] = Inf * ones(T, nx * ny) # first part equal to infinity
+    ucon[nx * ny + 1 : end] = zeros(T, 2 * nx + 2 * ny) # second part part equal to zero
+
+    v = vec(v_D) #convert to vector
+
+    nlp = ADNLPModel(Objective, v, constraints, lcon, ucon)
+    return nlp
+end 
+
+

src/ADNLPProblems/rocket.jl

@@ -0,0 +1,91 @@
+#  Find the surface with minimal area, given boundary conditions, 
+#  and above an obstacle.
+
+#  This is problem 10 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)
+#  classification OBR2-AN-V-V
+
+function rocket(; n::Int = default_nvar, type::Val{T} = Val(Float64), kwargs...) where {T}
+
+    # Initialisation
+    # Constants
+    h_0 = T(1) # height initialization   
+    v_0 = T(0) # speed initialization   
+    m_0 = T(1) # mass initialization  
+    g_0 = T(1) # gravity initialization
+
+    # Parameters
+
+    h_c = T(500) # Used for drag
+    v_c = T(620) # Used for drag
+    m_c = T(0.6) # Fraction of initial mass left at end
+
+    c = T(1/2 * (g_0*h_0)^2)        # Thrust-to-fuel mass
+    m_f = T(m_0 * m_c)              # final mass
+    T_max = T(3.5 * g_0 * m_0)      # maximal rocket thrust
+    D_c = T(1/2 * v_c * (m_0/g_0))  # Drag scaling
+
+    function f(X)        #Objective function
+        S = eltype(X)
+
+        v = zeros(S, n)  # velocity vector
+        h = zeros(S, n)  # height vector
+        g = zeros(S, n)  # gravity vector
+        m = zeros(S, n)  # mass vector
+        D = zeros(S, n)  # drag vector
+
+        v[1] = S(v_0)          # velocity vector initialization
+        h[1] = S(h_0)          # height vector initialization
+        g[1] = S(g_0)          # gravity vector initialization
+        m[1] = S(m_0)          # mass vector initialization
+        D[1] = S(D_c*(v_0^2))  # drag vector initialization
+        for k=2:n
+            m[k] = S(m[k - 1] - Δt * X[k - 1] / c)                                  # update mass vector   
+            v[k] = S(v[k - 1] + Δt *((X[k - 1] - D[k - 1]) / m[k - 1] - g[k - 1]))  # update speed vector
+            h[k] = S(h[k - 1] + Δt * v[k - 1])                                      # update height vector
+            D[k] = S(D_c*(v[k]^2)*exp(-h_c*(h[k]-h_0)/h_0))                         # update drag vector
+            g[k] = S(g_0*(h_0/h[k])^2)                                              # update gravity vector
+        end
+        opt = -h[end]
+        return opt
+
+    end
+
+    function constraints(X)
+
+        S = eltype(X)
+
+        v = zeros(S, n)
+        h = zeros(S, n)
+        g = zeros(S, n)
+        m = zeros(S, n)
+        D = zeros(S, n)
+
+        v[1] = S(v_0)          # velocity vector initialization
+        h[1] = S(h_0)          # height vector initialization
+        g[1] = S(g_0)          # gravity vector initialization
+        m[1] = S(m_0)          # mass vector initialization
+        D[1] = S(D_c*(v_0^2))  # drag vector initialization
+        for k=2:n
+            m[k] = S(m[k - 1] - Δt * X[k - 1] / c)
+            v[k] = S(v[k - 1] + Δt *((X[k - 1] - D[k - 1]) / m[k - 1] - g[k - 1]))
+            h[k] = S(h[k - 1] + Δt * v[k - 1])
+            D[k] = S(D_c*(v[k]^2)*exp(-h_c*(h[k]-h_0)/h_0))
+        end
+        constraints = vcat(v, h .- h[1], m .- m_f)  # constraint vector for velocity, height, mass, thrust
+        return constraints
+
+    end
+    Δt  = T(1/(n-1)) 
+    X0 = T_max/2 * ones(T, n) 
+    lvar = zeros(T, n) 
+    uvar =  T_max/2 * ones(T, n) 
+    lcon = zeros(T, 3 * n)
+    ucon = T[i ≤ 2n ? T(Inf) : ( T(m_0 - m_f)) for i=1:3n]
+
+    return ADNLPModel(f, X0, lvar, uvar, constraints, lcon, ucon)
+end
+
+

src/Meta/minimalsurface.jl

@@ -0,0 +1,26 @@
+minimalsurface_meta = Dict(
+    :nvar => 400,
+    :variable_nvar => false,
+    :ncon => 480,
+    :variable_ncon => false,
+    :minimize => true,
+    :name => "minimalsurface",
+    :has_equalities_only => false,
+    :has_inequalities_only => false,
+    :has_bounds => false,
+    :has_fixed_variables => false,
+    :objtype => :other,
+    :contype => :general,
+    :best_known_lower_bound => -Inf,
+    :best_known_upper_bound => Inf,
+    :is_feasible => missing,
+    :defined_everywhere => missing,
+    :origin => :unknown,
+
+)
+get_minimalsurface_nvar(; n::Integer = default_nvar, kwargs...) = 400
+get_minimalsurface_ncon(; n::Integer = default_nvar, kwargs...) = 480
+get_minimalsurface_nlin(; n::Integer = default_nvar, kwargs...) = 0
+get_minimalsurface_nnln(; n::Integer = default_nvar, kwargs...) = 480
+get_minimalsurface_nequ(; n::Integer = default_nvar, kwargs...) = 80
+get_minimalsurface_nineq(; n::Integer = default_nvar, kwargs...) = 400

src/Meta/rocket.jl

@@ -0,0 +1,25 @@
+rocket_meta = Dict(
+  :nvar => 100,
+  :variable_nvar => true,
+  :ncon => 300,
+  :variable_ncon => true,
+  :minimize => true,
+  :name => "rocket",
+  :has_equalities_only => false,
+  :has_inequalities_only => true,
+  :has_bounds => true,
+  :has_fixed_variables => false,
+  :objtype => :other,
+  :contype => :general,
+  :best_known_lower_bound => -Inf,
+  :best_known_upper_bound => Inf,
+  :is_feasible => missing,
+  :defined_everywhere => missing,
+  :origin => :unknown,
+)
+get_rocket_nvar(; n::Integer = default_nvar, kwargs...) = 1 * n + 0
+get_rocket_ncon(; n::Integer = default_nvar, kwargs...) = 3 * n + 0
+get_rocket_nlin(; n::Integer = default_nvar, kwargs...) = 0
+get_rocket_nnln(; n::Integer = default_nvar, kwargs...) = 3 * n + 0
+get_rocket_nequ(; n::Integer = default_nvar, kwargs...) = 0
+get_rocket_nineq(; n::Integer = default_nvar, kwargs...) = 3 * n + 0