Трудности достижения сходимости с математической оптимизацией NLOPT для физической проблемы ⇐ C#

[Anonymous]

Я кодирую физический симулятор для моделирования статического поведения объектов, плавающих в воде (в основном корабли). Положение объекта описывается глубиной погружения (draft) и два угла (heel и trim ). На объекте всегда есть как минимум 2 силы - вес и плавучесть, но может быть больше. Сила описывается как величина, происхождение и вектор. Симулятор должен сбалансировать силы, чтобы чистая сила вверх и поворотные моменты (для каблуки и обрезки) были почти озвучиваемыми. < /P>
в основном проблема математической оптимизации, поэтому я использую библиотеку NLOPT (в частности Nloptnet, так как это проект C#). Существует три параметра управления (черновик, каблук и отделка) и три цели, чтобы минимизировать (netUpwardForce, Heelmoment и Trimmoment ). Изменение значения одного параметра влияет на все три цели, но всегда больше, чем другие (например, изменение каблука влияет на каблук гораздо больше, чем другие цели). К сожалению, NLOPT не допускает несколько целей, поэтому их нужно нормализовать и объединить в один. < /P>
Вот код: < /p>
double OptimizeFunc(double[] parameters) //parameters: Draft, heel and trim
{
Model.SetPosition(parameters[0], parameters[1], parameters[2]);

List 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)
}

Проблема, которая у меня есть, заключается в том, что nloptsolver не может минимизировать цель (найдите черновик , каблук и trim значения, которые приводят к почти нулевая сеть вверх и моменты пятки/отделки). Для конкретного теста, который я работаю, это график из 50 попыток, которые делает Nloptsolver, построив оценку (значение, возвращаемое оптимизациейфунком ), а также проект , Heel и trim значения, которые он использовал. /p>
Несмотря на то, что решатель считает, что он нашел локальный минимум, путем вручную настраивая черновик , каблук и trim Значения Я знаю, что на самом деле лучшее решение - это Проект = 0,202 , heel = 0 , trim = -0,32 , что приводит к оценке = 0,001 (и это может быть усовершенствовано еще больше, но я остановился там) Полем Между тем, как вы можете видеть на графике, Solver только удается достичь оценки = 0,013 , что даже не близко. Я прошел дополнительные тесты, чтобы наметить все пространство оптимизации, и подтвердил, что нет других минимумов, кроме этого. Я также пытался ограничить границы очень близко вокруг этого известного минимума, но решатель все еще не может найти его. > очень чувствительны к NetUpwardforce не нулевым, поэтому я попытался добавить весовых коэффициентов, чтобы оптимизатор приоритет минимизации NetUpwardforce сначала, но это не помогает. Я также попытался отказаться от оптимизации для каблучения и Trummoment , а вместо этого установил оптимизатор, чтобы минимизировать горизонтальные расстояния между центров плавучести и сил веса, но это тоже не помогло. Я также попытался поднять NetUpforComp , HeelMomMomcomp и trimmomcomp к более высоким силам, и это немного помогло, но все еще не решает проблему.
< P> Может ли кто -нибудь испытать с математической оптимизацией предположить, что может быть не так с моим подходом? Это способ, которым я вычисляю Total в OptimizeFunc ? Или я должен использовать другой алгоритм вместо бобивы? Я попробовал несколько, но ни один не дал лучших результатов.draft = 0.202, Heel = 0 , trim = -0,32 ) в качестве начальных начальных значений для NLOPT, тогда он блокируется к этому решению и принимает его как действительное после нескольких циклов. Но это не так, когда начинается с нулевых значений, даже если они относительно очень близки к правильным значениям решения. Полем Я уверен, что они работают правильно. < /P>
public static class StabilityCalculator
{
public static StabilityCalculatorData CalculateStability(List 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
);
}
}


Подробнее здесь: https://stackoverflow.com/questions/793 ... n-for-a-ph