Optimising a non-linear user-defined function using NLopt

Dear Julia experts,

This is my second Julia code. I have a physical problem to solve in which a waving filament in fluid medium propels itself (think of a beating sperm tail).

I have a MATLAB code which optimises the shape-modes of a waving filament for the maximum propulsion speed in the form of an efficiency term. The 3 functions sgf_3d_fs, elementintegrate and threenestedloops are needed for the main function filamentefficiency. The filamentefficiency is a function of a vector A which contains the shape modes. So for a given vector A, it will return the efficiency (Effc) of the waving filament.

For example, for just a sine wave, the efficiency of the waving filament is

julia> Effc = filamentefficiency([1.0, 0.0, 0.0, 0.0, 0.0])

and for just a wave that looks like a saw-tooth, the efficiency becomes higher,

julia> Effc = filamentefficiency([0.9653, -0.1037, 0.0351, -0.0164, 0.0088])

I got the latter optimised-modes from a MATLAB program using fmincon (sequential quadratic programming) function. I am trying to solve the same problem in Julia because MATLAB is very slow when I want to optimise efficiency with modes > 15. I am trying to use NLopt.jl to maximise the function. I have tried using JuMP.jl as well with NLopt.jl but I could not make it work.

I get a strange output when I run the program below. I have tried derivative-free algorithms as well but it always stops after 1 iterations and gives either 0 or Inf. Clearly, I am not using NLopt correctly. I will be grateful if someone can help me figure this out. The last 5 lines of the code are relevant for the optimisation part using NLopt.

got -Inf at [1.0, 0.0, 0.0, 0.0, 0.0] after 1 iterations (returned FORCED_STOP)
using FastGaussQuadrature, LinearAlgebra, BenchmarkTools, TimerOutputs, StaticArrays, NLopt

function sgf_3d_fs(x,y,z,x0,y0,z0)
    dx = x-x0
    dy = y-y0
    dz = z-z0
    dxx = dx^2
    dxy = dx*dy
    dxz = dx*dz
    dyy = dy^2
    dyz = dy*dz
    dzz = dz^2

    r  = sqrt(dxx+dyy+dzz);
    r3 = r*r*r
    ri  = 1.0/r
    ri3 = 1.0/r3

    g11 = ri  + dxx*ri3
    g12 = dxy*ri3
    g13 = dxz*ri3

    g21 = g12
    g22 = ri + dyy*ri3
    g23 = dyz*ri3

    g31 = g13
    g32 = g23
    g33 = ri + dzz*ri3

    GEk = @SMatrix [ g11 g12 g13
                     g21 g22 g23
                     g31 g32 g33 ]
    return GEk

function elementintegrate(xc,yc,tcx,tcy,sc,x1,y1,x2,y2,s1,s2,xg,wg,ng)
    # x1,y1,x2,y2 are coordinates of element j with nodes at j and j+1
    # xc, yc are the coordinates of collocation point
    GE = zeros(SMatrix{3,3,Float64}) #
    GEself = zeros(SMatrix{3,3,Float64}) #
    d = 0.001 # regularisation parameter
    for k = 1:ng
        xi = xg[k]
        xint= 0.5*( (x2 + x1) + (x2 - x1)*xi)
        yint = 0.5*( (y2 + y1) + (y2 - y1)*xi)
        sint = 0.5*( (s2 + s1) + (s2 - s1)*xi)
        hl = 0.5*sqrt( (x2-x1)^2 + (y2-y1)^2)
       # Call Green's function for Stokes equation
        GEk = sgf_3d_fs(xint,yint,0,xc,yc,0) # two dimensional problem: z=0
        GE += GEk*hl*wg[k]
        ds = abs(sc-sint)
        dss = sqrt(ds^2 + d^2)
        Gs11 = ((1.0 + tcx*tcx)/dss)*hl*wg[k]
        Gs12 = ((0.0 + tcx*tcy)/dss)*hl*wg[k]
        Gs21 = ((0.0 + tcy*tcx)/dss)*hl*wg[k]
        Gs22 = ((1.0 + tcy*tcy)/dss)*hl*wg[k]
        GEself += @SMatrix [ Gs11 Gs12 0
                             Gs21 Gs22 0
                             0    0    0 ]
    return GE, GEself

function threenestedloops(xm,ym,tmx,tmy,sm,vm,ls,s,c,x,y,xg,wg,ng,nx)
   # Initialize the matrix and right hand side
    RM = zeros(2*(nx-1) + 1,2*(nx-1) + 1)
    rhs = zeros(2*(nx-1) + 1,1)
   #Threads.@threads for j = 1:nx-1   # use Threads if needed
    for j = 1:nx-1 # looping over columns first
        j1 = (j-1)*2 + 1
        j2 = (j-1)*2 + 2

        tcx = tmx[j]
        tcy = tmy[j]

        # Local operator terms( \Lambda [f_h] )
        RM[j1,j1] += -(1*(2-c) - tcx*tcx*(c+2))/8π
        RM[j1,j2] += -(0*(2-c) - tcx*tcy*(c+2))/8π
        RM[j2,j1] += -(0*(2-c) - tcy*tcx*(c+2))/8π
        RM[j2,j2] += -(1*(2-c) - tcy*tcy*(c+2))/8π

        rhs[j1] = 0
        rhs[j2] = vm[j]
        RM[2*(nx-1) + 1,j1] = ls[j]
        RM[j1,2*(nx-1) + 1] = 1

        for i = 1:nx-1
            i1 = (i-1)*2 + 1
            i2 = (i-1)*2 + 2
            xc = xm[i]
            yc = ym[i]
            tcx = tmx[i]
            tcy = tmy[i]
            sc = sm[i]
            # Single element k-Loop Integration using quadratures
            GE,GEself = elementintegrate(xc,yc,tcx,tcy,sc,x[j],y[j],x[j+1],y[j+1],s[j],s[j+1],xg,wg,ng)

            # Non-local operator terms( K [f_h] )

            RM[i1,j1] += GE[1,1]/8π
            RM[i1,j2] += GE[1,2]/8π
            RM[i2,j1] += GE[2,1]/8π
            RM[i2,j2] += GE[2,2]/8π

            RM[i1,i1] += -GEself[1,1]/8π
            RM[i1,i2] += -GEself[1,2]/8π
            RM[i2,i1] += -GEself[2,1]/8π
            RM[i2,i2] += -GEself[2,2]/8π
    return RM, rhs

function filamentefficiency(A)
    # Specify the number of nodes
    nx = 101
    # Note that the number of elements = nx-1
    # How many points needed for integration using Gauss Quadrature?
    ng = 20
    # Specify the odd-numbered mo20
    # A = [0.9653, -0.1037, 0.0351, -0.0164, 0.0088]
    # time
    t = 0
    # Points and weights "using FastGaussQuadrature"
    xg, wg = gausslegendre(ng)
    # Generate a sin curve
    # x and y are nodes while x_m and y_m are mid points of an element
    x = LinRange(-π,π,nx)
    nm = 1:1:length(A)
    y = (sin.((x.-t).*(2nm'.-1)))*A #(Only odd numbered modes)
    # Mid-points (collocation points) and tangents at the mid points.
    xm = 0.5*(x[2:nx] .+ x[1:nx-1])
    ym = 0.5*(y[2:nx] .+ y[1:nx-1])
    ls = sqrt.((x[2:nx] .- x[1:nx-1]).^2 .+ (y[2:nx] .- y[1:nx-1]).^2) # length of each element
    tmx = (x[2:nx] .- x[1:nx-1])./ls
    tmy = (y[2:nx] .- y[1:nx-1])./ls
    sqrt.(tmx.^2 + tmy.^2)
    # Arc length of node points
    s = zeros(nx,1)
    for i = 1:nx-1
        s[i+1] = s[i] + ls[i]
    # Arc length of mid points
    sm = 0.5*(s[2:nx] .+ s[1:nx-1])
   # Filament velocity at mid points (forcing for the equations)
    vm = (-(2nm'.-1).*cos.((xm.-t).*(2nm'.-1)))*A
    ρ = 0.001
    c = 2*log(ρ/s[nx]) + 1
    RM, rhs = threenestedloops(xm,ym,tmx,tmy,sm,vm,ls,s,c,x,y,xg,wg,ng,nx)
    fh = RM\rhs
    Effc = - (fh[2*(nx-1)+1])^2/(sum(fh[2:2:2*(nx-1)].*ls.*vm))
    return Effc
    # println("The x-velocity is ", fh[2*(nx-1)+1])
    #display(plot(x,y, aspect_ratio=:equal))
    #display(quiver!(xm,ym,quiver=(tmx,tmy), aspect_ratio=:equal))

modes = 5
opt = Opt(:LD_SLSQP, modes)
opt.xtol_rel = 1e-4
opt.max_objective = filamentefficiency
(maxf,maxx,ret) = optimize(opt,  [1.0, 0.0, 0.0, 0.0, 0.0])
numevals = opt.numevals # the number of function evaluations
println("got $maxf at $maxx after $numevals iterations (returned $ret)")

In NLopt, even if you don’t use grad argument your objective function is expected to get 2 arguments: an optimization variable and its gradient. However, you are not expected to fill grad when you use non-gradient based methods. I had the same issue several times and it hits hard.

1 Like

Thank you tomaklutfu. Sorry, I don’t understand what is the solution to this issue then. You are right that I don’t have a gradient available as in the tutorial of NLopt github page. I have used derivative-free methods as well but I still have the same problem. So how do I proceed ?

Make your function accept 2 arguments. Even if the method derivative-free nlopt uses 2 argument objective functions. In your function just ignore the second argument.

function filamentefficiency(A, grad) #Ignore the second argument if method is derivative-free
1 Like

Thank you man. You are a life-saver. I changed the function as you said. But opt = Opt(:LD_SLSQP, modes) does not work and quite understandably as it needs gradient. So I changed the algorithm to opt = Opt(:LN_COBYLA, modes) and I got the optimised modes.

got 0.011027008689477235 at [0.9574894546716676, -0.11139750102499951, 0.03920170352882533, -0.01879630947003879, 0.010659052201069571]

SQP works in MATLAB without me specifying any gradient, so it must be doing something to generate a gradient. I should try finding out what exactly it does.

On a different note, and taking a step back, what other packages I can try for this problem? Perhaps JuMP with a nonlinear optimisation algorithm? However, I cannot find any mention of derivative-free algorithms there. So, one has to make use of Automatic differentiation and also do something special for making it take vector inputs User-defined functions with vector inputs

For the fast prototyping and without making much changes, you can use FiniteDifference.jl for derivatives and test it. I do not know much about JuMP, I have not used it personally.

For AD to work, hard typed type information should go and it should depend on inputs element type. It is only this 2 allocations, I think, should only change.