Difficulties achieving convergence with NLopt mathematical optimization for a physics problem

Question

I am coding a physics simulator for simulating the static behavior of objects floating in water (mainly ships). Object position is described by immersion depth (draft), and two angles (heel and trim). There are always at least 2 forces acting on the object - weight and buoyancy, but there may be more. Force is described as a magnitude, origin and vector.The simulator has to balance the forces so that the net upwards force and turning moments (for heel and trim) are near-zero.

This is basically a mathematical optimization problem, so I'm using NLopt library (specifically NLoptNet since this is a C# project). There are three control parameters (draft, heel and trim), and three objectives to minimize (netUpwardForce, heelMoment and trimMoment). Changing value of one parameter affects all three objectives, but one always more than others (for example, changing heel affects heelMoment much more than other objectives). Unfortunately, NLopt does not allow for multiple objectives, so they have to be normalized and combined into one.

Here is the code:

double OptimizeFunc(double[] parameters) //parameters: Draft, heel and trim
{
    Model.SetPosition(parameters[0], parameters[1], parameters[2]);

    List<Force> forces = new(){ Model.GetBuoyantForce(), Model.GetWeightForce() };

    StabilityData stabilityData = StabilityCalculator.Calculate(forces); //calculates net upward force, heeling moment and trimming moment

    //Values for this specific test
    const double maxTheoreticalBuoyancy = 26000; 
    const double maxTheoreticalHeelingMoment = 61000;
    const double maxTheoreticalTrimmingMoment = 128000;

    double netUpForceComp = Math.Abs(stabilityData.NetUpwardForce / maxTheoreticalBuoyancy);
    double heelMomComp    = Math.Abs(stabilityData.HeelingMoment  / maxTheoreticalHeelingMoment);
    double trimMomComp    = Math.Abs(stabilityData.TrimmingMoment / maxTheoreticalTrimmingMoment));

    double total = netUpForComp + heelMomComp + trimMomComp;    

    return total;
}

var parameters = { 0, 0, 0 } ;

using (solver = new NLoptSolver(NLoptAlgorithm.LN_BOBYQA, (uint)parameters.Count))
    {
        solver.SetLowerBounds(new double[] {0, -180 , -90});
        solver.SetUpperBounds(new double[] {Model.GetModelHeight(), 180 , 90});
        solver.SetMinObjective(OptimizeFunc);
        
        var nLoptResult = solver.Optimize(new double[] {0, 0 , 0}, out var finalScore);
        
        Debug.Print(finalScore)
    }

The issue I'm having is that the NLoptSolver fails to minimize the objective (find the draft, heel and trim values that result in near-zero net upward force and heel/trim moments). For a specific test I'm running, this is the graph of 50 attempts that the NLoptSolver makes, plotting the score (value returned by OptimizeFunc), as well as draft, heel and trim values it used.

Even though the solver thinks it found a local minimum, by manually tweaking draft, heel and trim values I know that the actually best solution is draft = 0.202, heel = 0, trim = -0.32, which results in a score = 0.001 (and it could be refined even more, but I stopped there). Meanwhile, as you can see from the graph, solver only manages to achieve score = 0.013, which is not even close. I ran additional tests to map out the entire optimization space, and confirmed that there are no other minimums besides this one. I also tried to restrict bounds very close around this known minimum, but the solver still fails to find it.

I suspect that the issue is that HeelingMoment and TrimmingMoment are very sensitive to NetUpwardForce not being zero, so I tried adding weighting factors so that the optimizer prioritizes minimizing NetUpwardForce first, but that doesn't help matters. I also tried to discard optimizing for HeelingMoment and TrimmingMoment, and instead set the optimizer to minimize the horizontal distances between the buoyancy and weight force centers, but that didn't help either. I also tried to raise netUpForComp, heelMomComp and trimMomComp to higher powers, and that helped a little, but still doesn't solve the issue.

Can anyone experienced with mathematical optimization suggest what could be wrong with my approach? Is it the way I compute the total in OptimizeFunc? Or should I be using another algorithm instead of BOBYQA? I tried several, but none gave any better results.

EDIT: additional detail: if I pass the correct solution parameters (draft = 0.202, heel = 0, trim = -0.32) as the initial starting values for Nlopt, then it does lock on to that solution and accepts it as valid after a few cycles. But it doesn't when starting from zero values, even though they are relatively very close to the correct solution values.

EDIT 2: as per request, I am adding the code of supporting classes. I'm pretty sure these work correctly though.

public static class StabilityCalculator
{
    public static StabilityCalculatorData CalculateStability(List<Force> forces)
    {        
        Coordinate center = forces[0].Origin
        
        Force netUpwardForce = forces.Sum(f => f.Magnitude * f.Vector.Y);

        Coordinate totalMoment = new(0, 0, 0);
        foreach (Force force in forces)
        {
            Coordinate relativePosition = force.Origin - center;
            Coordinate momentVector = Coordinate.CrossProduct(relativePosition, force.Vector * force.Magnitude);
            totalMoment += momentVector;
        }

        return new StabilityCalculatorData
        {
            HeelingMoment = totalMoment.Z, 
            YawingMoment = totalMoment.Y,  
            TrimmingMoment = totalMoment.X,
            NetUpwardForce = netUpwardForce
        };
    }
}

public class StabilityCalculatorData
{
    public double HeelingMoment { get; set; } // Moment about X-axis
    public double TrimmingMoment { get; set; } // Moment about Z-axis
    public double YawingMoment { get; set; } // Moment about Y-axis
    public double NetUpwardForce { get; set; } // Net upward force in Newtons
}

public class Force(double magnitude, Coordinate origin, Coordinate vector)
{
    public double Magnitude { get; } = magnitude;
    public Coordinate Origin { get;  } = origin;
    public Coordinate Vector { get; } = vector.Normalize();
    
    public static Force? Combine(Force? f1, Force? f2)
    {
        if (f1 == null && f2 != null) { return f2;}
        if (f2 == null && f1 != null) { return f1;}
        if (f1 == null && f2 == null) { return null; }

        var vector1 = new Coordinate(f1!.Vector.X * f1.Magnitude, f1.Vector.Y * f1.Magnitude, f1.Vector.Z * f1.Magnitude);
        var vector2 = new Coordinate(f2!.Vector.X * f2.Magnitude, f2.Vector.Y * f2.Magnitude, f2.Vector.Z * f2.Magnitude);

        var resultantVector = vector1 + vector2;
        var resultantMagnitude = resultantVector.Magnitude();

        var normalizedResultantVector = resultantVector.Normalize();

        return new Force(resultantMagnitude, f1.Origin, normalizedResultantVector);
    }
}

public class Coordinate(double x, double y, double z)
{
    public double X { get; set; } = x;
    public double Y { get; set; } = y;
    public double Z { get; set; } = z;

    public static Coordinate operator +(Coordinate c1, Coordinate c2)
    {
        return new Coordinate(c1.X + c2.X, c1.Y + c2.Y, c1.Z + c2.Z);
    }

    public static Coordinate operator -(Coordinate c1, Coordinate c2)
    {
        return new Coordinate(c1.X - c2.X, c1.Y - c2.Y, c1.Z - c2.Z);
    }

    public static Coordinate operator *(Coordinate c, double scalar)
    {
        return new Coordinate(c.X * scalar, c.Y * scalar, c.Z * scalar);
    }

    public static Coordinate operator /(Coordinate c, double scalar)
    {
        return new Coordinate(c.X / scalar, c.Y / scalar, c.Z / scalar);
    }
    
    public Coordinate Normalize()
    {
        var magnitude = Math.Sqrt(X * X + Y * Y + Z * Z);
        return new Coordinate(X / magnitude, Y / magnitude, Z / magnitude);
    }

    public double Magnitude()
    {
        return Math.Sqrt(X * X + Y * Y + Z * Z);
    }

    public static Coordinate CrossProduct(Coordinate c1, Coordinate c2)
    {
        return new Coordinate(
            c1.Y * c2.Z - c1.Z * c2.Y,
            c1.Z * c2.X - c1.X * c2.Z,
            c1.X * c2.Y - c1.Y * c2.X
        );
    }
}

Answer 1

Being derivative free, this solver must estimate an appropriate step size. The solver has bounds of order 360, and a 0 initial guess, so it may be having great difficulty finding an appropriate step size heuristically.

You could try providing a typical step size as described in the docs. https://nlopt.readthedocs.io/en/latest/NLopt_Reference/#initial-step-size.

So it can take steps which have the resolution to detect the solution.

(A better initial guess will always help also).