Мобильная версияЭто по сути, это задача математической оптимизации, поэтому я использую библиотеку NLopt (в частности, NLoptNet, поскольку это проект C#). Существует три параметра управления (осадка, крен и дифферент) и три цели, которые необходимо минимизировать (netUpwardForce, heelMoment и TrimMoment). Изменение значения одного параметра влияет на все три цели, но одна всегда больше, чем другие (например, изменение heel влияет на heelMoment гораздо сильнее, чем другие цели). К сожалению, NLopt не позволяет использовать несколько целей, поэтому их приходится нормализовать и объединить в одну.
Вот код:
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 не может минимизировать цель (найти значения уклона, крена и обрезки, которые приводят к почти нулевая чистая восходящая сила и моменты крена/дифферента). Для конкретного теста, который я провожу, это график 50 попыток, который делает NLoptSolver, на котором отображается оценка (значение, возвращаемое OptimizeFunc), а также уклон, пятка и обрезать использованные значения.
Несмотря на то, что решатель думает, что нашел локальный минимум, путем ручной настройки значений уклона, крена и обрезки я знаю, что на самом деле лучшее решение – черновик = 0,202, пятка = 0, обрезка = -0,32, что дает оценку = 0,001 (и это можно было бы еще уточнить, но я остановился на этом). Между тем, как видно из графика, решателю удается достичь только оценки = 0,013, что даже близко не близко. Я провел дополнительные тесты, чтобы отобразить все пространство оптимизации, и подтвердил, что других минимумов, кроме этого, нет. Я также пытался ограничить границы очень близко к этому известному минимуму, но решатель все равно не смог его найти.
Я подозреваю, что проблема в том, что HeelingMoment и TrimmingMoment очень чувствительны к тому, что NetUpwardForce не равен нулю, поэтому я попробовал добавить весовые коэффициенты, чтобы оптимизатор сначала отдавал приоритет минимизации NetUpwardForce, но это не помогает. Я также пытался отказаться от оптимизации HeelingMoment и TrimmingMoment и вместо этого настроить оптимизатор на минимизацию горизонтальных расстояний между центрами сил плавучести и веса, но это тоже не помогло. Я также пытался повысить степень netUpForComp, heelMomComp и TriMomComp, и это немного помогло, но проблему все равно не решило.
Может ли кто-нибудь, имеющий опыт математической оптимизации, подсказать, что может быть не так с моим подходом? Это способ вычисления суммы в OptimizeFunc? Или мне следует использовать другой алгоритм вместо BOBYQA? Я попробовал несколько, но ни один из них не дал лучших результатов.
РЕДАКТИРОВАТЬ: дополнительная информация: если я передам правильные параметры решения (draft = 0.202, heel = 0, Trim = -0,32) в качестве начальных начальных значений для Nlopt, затем он фиксируется на этом решении и принимает его как допустимое после нескольких циклов. Но этого не происходит при запуске с нулевых значений, хотя они относительно очень близки к правильным значениям решения.
РЕДАКТИРОВАТЬ 2: по запросу я добавляю код вспомогательных классов. . Хотя я почти уверен, что они работают правильно.
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
