Solve FORCED STOP issues

[I-Fons]

Dear all,

I would appreciate some help with this problem: I need to optimize the cost_fcn under a constraint on the yield loss (see yield_loss_v2), which must be between 0 and 0.95. When I run the code I receive a FORCED_STOP ret code and I don't know why.

module CaseStudy0

# add DiffEqFlux DiffEqOperators DiffEqSensitivity DifferentialEquations GalacticOptim NumericalIntegration Random Distributions

using DiffEqFlux
using DiffEqOperators
using DiffEqSensitivity
using DifferentialEquations
using GalacticOptim
using BlackBoxOptim
using NumericalIntegration
using Random
using Distributions
using MathOptInterface
using NLopt

# include("test_ICS.jl")

greet() = print("Hello World!")

mutable struct column_parameters
    Ac 
    ϵ 
    φ
    dp
    L
end

mutable struct LKM_parameters
    km_max 
    S1 
    S2 
    qsat
    H
end

mutable struct parameter_wrapper
    weights
    Δx
    n
    pLKM
    pCol
    c_in
    u0
    p
end

function km(km_max,S1,S2,qsat, q )
    km_max*(S1 + (1.0 - S1)*(1.0-q/qsat)^S2)
end

function q_eq(H,qsat, c ) 
    (c*H)/(1.0+c*H/qsat)
end

function TBC(pCol , c_in , pLKM)
    # Evaluation of the total binding capacity of the column
    # DOUBT: Use total column volume or only "liquid" volume
    (pCol.L*pCol.Ac*c_in*pLKM.H)/(1+c_in*pLKM.H/pLKM.qsat) 
end

function loaded_mass(optVar , parameters_all , M)
    # Make a prediction to compute loaded mass
    times = collect(0:optVar[2])
    cq = CaseStudy0.predict(parameters_all.Δx ,parameters_all.n ,optVar[1] ,parameters_all.c_in, parameters_all.pLKM , parameters_all.pCol , parameters_all.u0 , times ,  parameters_all.p)
    q_L_tend = cq[(parameters_all.n+1):end, end]
    # Trapezoidal rule used for numerical integration (further improvement)
    for i=2:parameters_all.n
        M = M + parameters_all.Δx * (q_L_tend[i-1]+q_L_tend[i])/2
    end
    return M
end

function yield_loss_v2(optVar, parameters_all)
    M = loaded_mass(optVar,parameters_all,0)
    capacity = TBC(parameters_all.pCol , parameters_all.c_in , parameters_all.pLKM)
    M/capacity - 0.95 # yield_loss = M/TBC <= 0.95
end

function cost_fcn(optVar , parameters_all)
    # optVar[1] is Q (feed flow rate) and optVar[2] is t_load (loading time)
    a_c = @view parameters_all.weights[1]
    a_m = @view parameters_all.weights[2]
    a_t = @view parameters_all.weights[3]
  
    M = loaded_mass(optVar,parameters_all,0)

    # Weighted sum of the three cost drivers
    a_c*parameters_all.pCol.L*parameters_all.pCol.Ac + a_m*M + a_t*optVar[2]
end

function PDE(du ,u ,p ,t ,Δx ,n ,Q ,c_in, pLKM , pCol)
    c = @view u[1:n]
    q = @view u[n+1:n*2]
    
    dc = @view du[1:n]
    dq = @view du[n+1:n*2]

    v = Q/pCol.Ac/pCol.ϵ
    D = pCol.dp*pCol.Ac*v/2 # Dispersion coefficient computed using reduced Van-Deemter-equation
    
    DB = CenteredDifference(2,2,Δx,n,D)
    VB = CenteredDifference(1,2,Δx,n,-v)
    bcc = RobinBC((1.0, 0.0, c_in), (0.0,1.0,0.0),Δx)

    @. dq = km(pLKM.km_max ,pLKM.S1 ,pLKM.S2 ,pLKM.qsat , q)*(q_eq(pLKM.H, pLKM.qsat, c) -  q)
    dc .= (DB+VB)*bcc*c .-pCol.φ*dq
end

function predict(Δx ,n ,Q ,c_in, pLKM , pCol , u0 , times ,  p)
    pde_closure = ODEFunction((du,u,p,t) -> PDE(du ,u ,p ,t ,Δx ,n ,Q ,c_in, pLKM , pCol))
    prob = ODEProblem(pde_closure,u0,(times[1],times[end]),p)
    sol = solve(prob,TRBDF2(),p=p,abstol = 1e-7, reltol = 1e-7,saveat= times,sensealg = ForwardDiffSensitivity())
    return Array(sol)
end

cb = function(p,l)
    println(l)
    if l < 1e-2
            return true
    else
            return false
    end
end

end # module

## Test the functions inside the module
# Column specifications
L = 10
dc = 0.5
Ac = pi*(dc/2)^2
V = L*pi*(dc/2)^2

ϵ = 0.36
dp =  8.5e-4
φ = (1-ϵ)/ϵ

pCol = CaseStudy0.column_parameters(Ac, ϵ , φ , dp , L)

# Parameters of the RHS of the constraints
b_yl = 0.95
b_r  = 0.95

# Parameters of the LKM
qsat = 83.8878 
H = 495.1227
S1 = 0.5169
S2 = 6.1372
km_max = 0.0541
# D = 1e-7

pLKM = CaseStudy0.LKM_parameters(km_max , S1 , S2 , qsat, H)

# Parameters of the stencil
n = 30 #spacepoints
Δx = L/(n-1)

# Guess values of the manipulated variables to initialize optimization
Q0 = 5
t_load0 = 250
optVar0 = [Q0 , t_load0]

# Other parameters to define
u0 = zeros(n*2) # Initial solution for the PDE system integrator
weights = [1,1,1] # weights associated to each cost driver
c_in = 0.5
p = "useless"
params_all = CaseStudy0.parameter_wrapper( weights , Δx ,n , pLKM , pCol , c_in , u0 , p)

## CONSTRAINED SOLUTION
# With NLopt
opt = CaseStudy0.Opt(:LN_COBYLA, 2)
opt.lower_bounds = [0., 0.]
opt.upper_bounds = [100, 1000]
opt.xtol_rel = 1e-4

opt.min_objective = CaseStudy0.cost_fcn
CaseStudy0.inequality_constraint!(opt, CaseStudy0.yield_loss_v2, 1e-3)

(minf,minx,ret) = CaseStudy0.NLopt.optimize(opt, optVar0)

Even when I try to solve the unconstrained problem I have the same FORCED_STOP...

Thanks in advance for the support

[I-Fons]

Duplicate of #127