# 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])
0.009672105882236201
``````

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])
0.011011146021855723
``````

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
end

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 ]
end
return GE, GEself
end

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π
end
end
return RM, rhs
end

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]
end
# 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))
end

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.

2 Likes