現在、 IPOPTとcsipoptを使用して C# で最適化の問題を実装しようとしています。
同じプロセスで同じデータを使用して同じ問題を何度も実行すると、毎回異なる結果になることがわかりました。とりわけ、CPU 時間はログ ファイルで増加し続け、最初の 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).
2 回目の反復ログ:
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 の専門家ではありません。
編集DisposeIPOPT インターフェイス
のメソッドは次のとおりです。
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;
}
}
これは問題自体を解放するように見えますが、ソルバーは解放しません。