Невоспроизводимые результаты с использованием IPOPT в C#
В настоящее время я пытаюсь реализовать проблему оптимизации в C#, используя IPOPT и csipopt.
Я заметил, что запуск одной и той же проблемы с одними и теми же данными снова и снова в одном и том же процессе каждый раз приводит к разным результатам. Среди прочего, время процессора увеличивается в файлах журнала, а начальный заголовок IPOPT печатается только в первом файле журнала. Однако, если я остановлю и перезапущу процесс, результаты будут воспроизводимы (т.е. итерация 0 одинакова каждый раз, итерация 1 тоже и т. Д.).
Журналы первой итерации:
List of options:
Name Value # times used
acceptable_constr_viol_tol = 0.01 1
acceptable_dual_inf_tol = 1e+010 1
acceptable_iter = 15 1
acceptable_obj_change_tol = 1e+020 1
acceptable_tol = 1e-006 1
alpha_for_y = primal 1
bound_mult_init_method = constant 1
bound_mult_init_val = 1 1
compl_inf_tol = 0.0001 2
constr_viol_tol = 0.0001 2
dual_inf_tol = 1 1
hessian_approximation = limited-memory 3
max_iter = 50 1
mu_init = 1 1
mu_strategy = monotone 2
mumps_scaling = 77 3
nlp_scaling_max_gradient = 100 1
recalc_y = no 1
recalc_y_feas_tol = 1e-006 0
tol = 1e-008 2
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
NOTE: You are using Ipopt by default with the MUMPS linear solver.
Other linear solvers might be more efficient (see Ipopt documentation).
This is Ipopt version 3.11.0, running with linear solver mumps.
Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 10
Number of nonzeros in Lagrangian Hessian.............: 0
Hessian approximation will be done in the space of all 9 x variables.
Scaling parameter for objective function = 1.000000e+000
objective scaling factor = 1
No x scaling provided
No c scaling provided
No d scaling provided
DenseVector "original x_L unscaled" with 9 elements:
[...]
Calling MUMPS-1 for symbolic factorization at cpu time 0.012 (wall 0.003).
Done with MUMPS-1 for symbolic factorization at cpu time 0.012 (wall 0.003).
MUMPS used permuting_scaling 5 and pivot_order 5.
scaling will be 77.
Calling MUMPS-2 for numerical factorization at cpu time 0.012 (wall 0.003).
Done with MUMPS-2 for numerical factorization at cpu time 0.012 (wall 0.003).
Number of doubles for MUMPS to hold factorization (INFO(9)) = 31
Number of integers for MUMPS to hold factorization (INFO(10)) = 164
Calling MUMPS-3 for solve at cpu time 0.012 (wall 0.003).
Done with MUMPS-3 for solve at cpu time 0.012 (wall 0.003).
Журналы второй итерации:
List of options:
Name Value # times used
acceptable_constr_viol_tol = 0.01 1
acceptable_dual_inf_tol = 1e+010 1
acceptable_iter = 15 1
acceptable_obj_change_tol = 1e+020 1
acceptable_tol = 1e-006 1
alpha_for_y = primal 1
bound_mult_init_method = constant 1
bound_mult_init_val = 1 1
compl_inf_tol = 0.0001 2
constr_viol_tol = 0.0001 2
dual_inf_tol = 1 1
hessian_approximation = limited-memory 3
max_iter = 50 1
mu_init = 1 1
mu_strategy = monotone 2
mumps_scaling = 77 3
nlp_scaling_max_gradient = 100 1
recalc_y = no 1
recalc_y_feas_tol = 1e-006 0
tol = 1e-008 2
This is Ipopt version 3.11.0, running with linear solver mumps.
Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 10
Number of nonzeros in Lagrangian Hessian.............: 0
Hessian approximation will be done in the space of all 9 x variables.
Scaling parameter for objective function = 1.000000e+000
objective scaling factor = 1
No x scaling provided
No c scaling provided
No d scaling provided
DenseVector "original x_L unscaled" with 9 elements:
[...]
Calling MUMPS-1 for symbolic factorization at cpu time 0.289 (wall 0.279).
Done with MUMPS-1 for symbolic factorization at cpu time 0.289 (wall 0.279).
MUMPS used permuting_scaling 5 and pivot_order 5.
scaling will be 77.
Calling MUMPS-2 for numerical factorization at cpu time 0.289 (wall 0.279).
Done with MUMPS-2 for numerical factorization at cpu time 0.289 (wall 0.279).
Number of doubles for MUMPS to hold factorization (INFO(9)) = 31
Number of integers for MUMPS to hold factorization (INFO(10)) = 164
Calling MUMPS-3 for solve at cpu time 0.289 (wall 0.279).
Done with MUMPS-3 for solve at cpu time 0.289 (wall 0.279).
Я предполагаю, что какой-то ресурс не освобождается при устранении моей проблемы, а некоторые значения хранятся где-то в памяти. Это становится большей проблемой, так как моя проблема иногда сходится, иногда не с теми же данными.
Моя проблема построена следующим образом:
using System;
using System.Linq;
using Cureos.Numerics;
namespace Ipopt
{
public class IpoptReconciliationProblem : IpoptProblem
{
private Double[] _standardDeviations;
private Double[] _averageValues;
private readonly JacobianHelper _jacobianHelper;
private readonly Func<Double[], Double[]> _constraints;
private Jacobian _jacobian;
private readonly Action<IterationResult> _callback;
private Boolean _refreshJac;
private Int32 _jacobianEvaluationFrequency;
public SolverResult Solve(Double[] averageValues, Double[] standardDeviations)
{
return this.Solve(averageValues, standardDeviations, averageValues.Clone() as Double[]);
}
public SolverResult Solve(Double[] averageValues, Double[] standardDeviations, Double[] initialValues)
{
Double obj;
this._averageValues = averageValues;
this._standardDeviations = standardDeviations;
var code = this.SolveProblem(initialValues, out obj);
var result = new SolverResult(initialValues, (Int32)code);
return result;
}
#region IpoptProblem methods override
public override IpoptBoolType eval_f(Int32 n, Double[] x, IpoptBoolType new_x, out Double obj_value,
IntPtr p_user_data)
{
//Objective function
obj_value = x.Select((d, i) => Math.Pow((d - this._averageValues[i]) / this._standardDeviations[i], 2)).Sum();
return IpoptBoolType.True;
}
public override IpoptBoolType eval_g(Int32 n, Double[] x, IpoptBoolType new_x, Int32 m, Double[] g,
IntPtr p_user_data)
{
//Contstraints
var updatedConstraints = this._constraints(x);
for (var i = 0; i < updatedConstraints.Length; i++)
{
g[i] = updatedConstraints[i];
}
return IpoptBoolType.True;
}
public override IpoptBoolType eval_grad_f(Int32 n, Double[] x, IpoptBoolType new_x, Double[] grad_f,
IntPtr p_user_data)
{
for (var i = 0; i < x.Length; i++)
{
grad_f[i] = 2 * ((x[i] - this._averageValues[i]) / this._standardDeviations[i]) *
(1 / this._standardDeviations[i]);
}
return IpoptBoolType.True;
}
public override IpoptBoolType eval_jac_g(Int32 n, Double[] x, IpoptBoolType new_x, Int32 m, Int32 nele_jac,
Int32[] iRow, Int32[] jCol, Double[] values, IntPtr p_user_data)
{
if (values == null)
{
for (var i = 0; i < iRow.Length; i++)
{
iRow[i] = this._jacobian.iRow[i];
}
for (var i = 0; i < jCol.Length; i++)
{
jCol[i] = this._jacobian.jCol[i];
}
}
else
{
if (this._refreshJac)
{
this._jacobian = this._jacobianHelper.UpdateJacobian(x, this._constraints, this._jacobian);
}
//else
//{
for (var i = 0; i < values.Length; i++)
{
values[i] = this._jacobian.values[i];
}
//}
}
return IpoptBoolType.True;
}
public override IpoptBoolType eval_h(Int32 n, Double[] x, IpoptBoolType new_x, Double obj_factor, Int32 m,
Double[] lambda, IpoptBoolType new_lambda, Int32 nele_hess, Int32[] iRow, Int32[] jCol, Double[] values,
IntPtr p_user_data)
{
return IpoptBoolType.True;
}
/// <summary>
/// Intermediate callback for each iteration.
/// </summary>
/// <param name="alg_mod">The algorithm mode.</param>
/// <param name="iter_count">The current iteration count.</param>
/// <param name="obj_value">The unscaled objective value at the current point.</param>
/// <param name="inf_pr">The unscaled constraint violation at the current point.</param>
/// <param name="inf_du">The scaled dual infeasibility at the current point.</param>
/// <param name="mu">log_10 of the value of the barrier parameter mu</param>
/// <param name="d_norm">The infinity norm (max) of the primal step.</param>
/// <param name="regularization_size">log_10 of the value of the regularization term for the Hessian of the Lagrangian in the augmented system. A dash (-) indicates that no regularization was done.</param>
/// <param name="alpha_du">The step size for the dual variables for the box constraints in the equality constrained formulation.</param>
/// <param name="alpha_pr">The step size for the primal variables x and lambda in the equality constrained formulation.</param>
/// <param name="ls_trials">The number of backtracking line search steps (does not include second-order correction steps).</param>
/// <param name="p_user_data">User data.</param>
/// <returns></returns>
public override IpoptBoolType intermediate(IpoptAlgorithmMode alg_mod, Int32 iter_count, Double obj_value,
Double inf_pr, Double inf_du,
Double mu, Double d_norm, Double regularization_size, Double alpha_du, Double alpha_pr, Int32 ls_trials,
IntPtr p_user_data)
{
this._callback(new IterationResult
{
ObjectiveFunction = obj_value,
IterationNumber = iter_count,
BacktrackingLines = ls_trials,
ConstraintViolation = inf_pr,
DualInfeasibility = inf_du,
DualVariablesStepSize = alpha_du,
Mu = mu,
PrimalStepNorm = d_norm,
PrimalVariablesStepSize = alpha_pr,
RegularizationTerm = regularization_size
});
this._refreshJac = iter_count > 0 && iter_count % this._jacobianEvaluationFrequency == 0;
return IpoptBoolType.True;
}
#endregion
#region Constructors
public IpoptReconciliationProblem(JacobianHelper jacobianHelper, Double[] lowerBounds,
Double[] upperBounds, Double[] constraintLowerBounds, Double[] constraintUpperBounds,
Func<Double[], Double[]> constraints, Jacobian jacobian, IpoptParams ipoptParams,
Action<IterationResult> callback)
: base(
lowerBounds.Length, lowerBounds, upperBounds, constraintLowerBounds.Length, constraintLowerBounds,
constraintUpperBounds, jacobian.iRow.Length, 0, true, true, true)
{
this._jacobianHelper = jacobianHelper;
this._constraints = constraints;
this._jacobian = jacobian;
this._callback = callback;
this.SetOptions(ipoptParams);
}
#endregion
private void SetOptions(IpoptParams ipoptParams)
{
this.AddOption("tol", ipoptParams.Tolerance);
this.AddOption("max_iter", ipoptParams.MaximumIterations);
this.AddOption("dual_inf_tol", ipoptParams.DualInfeasibilityTolerance);
this.AddOption("constr_viol_tol", ipoptParams.ConstraintViolationTolerance);
this.AddOption("compl_inf_tol", ipoptParams.ComplementaryInfTolerance);
// Acceptable solution
this.AddOption("acceptable_iter", ipoptParams.AcceptableIterationNumber);
this.AddOption("acceptable_tol", ipoptParams.AcceptableTolerance);
this.AddOption("acceptable_obj_change_tol", ipoptParams.AcceptableObjFctChange);
this.AddOption("acceptable_dual_inf_tol", ipoptParams.AcceptableDualInfeasibilityTolerance);
this.AddOption("acceptable_constr_viol_tol", ipoptParams.AcceptableConstraintViolationTolerance);
// Initialization
this.AddOption("bound_mult_init_val", ipoptParams.BoundMultiplier);
this.AddOption("bound_mult_init_method", ipoptParams.BoundMultiplierInitializationMethod);
this.AddOption("mu_init", ipoptParams.MuInit);
this.AddOption("mu_strategy", ipoptParams.MuStrategy);
//this.AddOption("nlp_upper_bound_inf", 1e-8);
//this.AddOption("bound_relax_factor", 1e-4);
//this.AddOption("nlp_scaling_min_value", 0.1);
//this.AddOption("bound_push", 1e-16);
//this.AddOption("bound_frac", 1e-16);
this.AddOption("alpha_for_y", ipoptParams.AlphaForY);
//this.AddOption("quality_function_max_section_steps",3);
this.AddOption("recalc_y_feas_tol", ipoptParams.RecalcYFeasTol);
this.AddOption("recalc_y", ipoptParams.RecalcY);
if (!String.IsNullOrEmpty(ipoptParams.LogFilePath))
{
this.OpenOutputFile(ipoptParams.LogFilePath, ipoptParams.LogLevel);
}
//JsonFile.Save(ipoptParams, @"d:\tmp\ipoptParams.json");
this.AddOption("nlp_scaling_max_gradient", ipoptParams.MaxGradient);
this.AddOption("mumps_scaling", ipoptParams.MumpsScaling);
this._jacobianEvaluationFrequency = ipoptParams.JacFreq;
}
}
}
Затем я называю это так (из консольного приложения):
var dic = results.Select((d, i) => new {Name = d.TagName, Index = i})
.ToDictionary(d => d.Name, d => d.Index);
var avg = results.Select(d => d.AverageValue).ToArray();
var stdDev = results.Select(d => d.StdDev).ToArray();
var lowerBounds = results.Select(d => d.Min).ToArray();
var upperBounds = results.Select(d => d.Max).ToArray();
var cLowerBounds = Enumerable.Range(0, constraints).Select(d => -Double.Epsilon).ToArray();
var cUpperBoundsBounds = cLowerBounds.Select(d => Double.Epsilon).ToArray();
var jacobianHelper = new JacobianHelper();
Func<Double[], Double[]> constraintCallback = values =>
{
var gasHpToComp = values[dic["Gas Outlet HP to Comp (kSm3/d)"]];
var gasLpToComp = values[dic["Gas Outlet LP to Comp (kSm3/d)"]];
var flareRecovery = values[dic["Flares Recovery"]];
var flareLP2 = dic.ContainsKey("Gas Comp LP2 to HP Flare (kSm3/d)")
? values[dic["Gas Comp LP2 to HP Flare (kSm3/d)"]]
: 0;
var gasBeforeTeg = gasHpToComp + gasLpToComp + flareRecovery + flareLP2;
var teg = values[dic["Gas TEG"]];
var r1 = gasBeforeTeg - teg;
var gasExport = values[dic["Gas Export"]];
var hp1Flare = values[dic["Gas Comp HP1 to HP Flare"]];
var fuelgas = values[dic["Fuel Gas"]];
var gasLift = values[dic["Gas Lift"]];
var gasAfterTeg = gasExport + hp1Flare + fuelgas + gasLift;
var r2 = teg - gasAfterTeg;
return new[] {r2, r1};
};
var jac = jacobianHelper.ComputeJacobian(avg, constraintCallback);
const String fileName = "ipoptLog";
var ipoptParams = new IpoptParams
{
LogLevel = 12
};
for (var i = 0; i < 5; i++)
{
ipoptParams.LogFilePath = Path.Combine(baseDirectory, $"{fileName}.{i}.txt");
SolverResult res;
using (var sut = new IpoptReconciliationProblem(jacobianHelper, lowerBounds, upperBounds, cLowerBounds,
cUpperBoundsBounds, constraintCallback, jac, ipoptParams, ir => { }))
{
res = sut.Solve(avg, stdDev);
}
Console.WriteLine("Resolution " + i + " - result: " + res.ReturnCode);
}
С другой стороны (но может быть связано?), Та же проблема сходится в модульном тесте и в консольном приложении, но не в службе Windows или в WebAPI.
Полный исполняемый код доступен в репозитории (VS2017), демонстрируя невоспроизводимые результаты.
Любая подсказка или руководство будет оценено. Я ни в коем случае не эксперт в вопросах оптимизации, тем более в IPOPT.
РедактироватьDispose метод в интерфейсе IPOPT выглядит следующим образом:
protected virtual void Dispose(bool disposing)
{
if (!m_disposed)
{
if (m_problem != IntPtr.Zero)
IpoptAdapter.FreeIpoptProblem(m_problem);
if (disposing)
{
m_problem = IntPtr.Zero;
}
m_disposed = true;
}
}
Это, кажется, освобождает саму проблему, но не решает.