NLopt Fails to Converge on Expected Minimum in C# Ship Stability Simulator

I'm developing a physics simulator to model the static behavior of ships floating in water. The simulation calculates the ship's position using three parameters: draft, heel, and trim. These affect the net upward force and the turning moments (heel and trim) generated by forces such as buoyancy and weight. The objective is to balance these forces so that the net upward force and moments are close to zero.

This is formulated as a mathematical optimization problem. I use the NLopt library (via NLoptNet in C#) with three normalization components combined into a single objective function. Changing one parameter affects all objectives, though each parameter has a dominant effect on one of them. Here’s the primary code snippet:

CSHARPCopied
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 problem is that while NLopt finds a local minimum (with a score around 0.013), manually testing suggests that an optimal solution exists at draft = 0.202, heel = 0, trim = -0.32 with a much lower score (approximately 0.001). Even restricting the parameter bounds around this known solution doesn’t help when starting from {0, 0, 0}.

I’ve tried several fixes including:
• Adjusting normalization weights and increasing the power of each term
• Prioritizing net upward force minimization separately
• Changing the objective to minimize horizontal distances between force centers

None of these approaches led NLopt to consistently find the narrow optimum unless I supplied the correct starting values. I suspect there may be an issue with the sensitivity of the combined objective or the algorithm’s ability to adjust step sizes given the large parameter bounds.

Has anyone experienced similar challenges with NLopt or have suggestions on adjusting the formulation, scaling, or algorithm choice (instead of BOBYQA) to reliably converge on the true optimum?

Additional supporting code for the force and stability calculations is included below:

CSHARPCopied
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
        );
    }
}