¿Cómo puedo emular la funcionalidad de Solver de Microsoft Excel (GRG no lineal) en C#?

Question

Tengo un problema de optimización no lineal con restricciones. Se puede resolver en Microsoft Excel con el complemento Solver, pero tengo problemas para replicarlo en C#.

Mi problema se muestra en el following spreadsheet. Estoy resolviendo el problema clásico A x = b pero con la advertencia de que todos los componentes de x deben ser no negativos. Entonces, en lugar de usar álgebra lineal estándar, utilizo Solver con la restricción no negativa, minimizando la suma de las diferencias al cuadrado, y obtengo una solución razonable. Intenté replicar esto en C# usando Microsoft Solver Foundation o Solver SDK. Sin embargo, parece que no puedo llegar a ninguna parte con ellos porque con MSF no puedo entender cómo definir el objetivo y con Solver SDK siempre obtengo el estado "óptimo" y una solución de todos los 0, que definitivamente no es ni siquiera local. mínimo.

Aquí está mi código de Solver SDK:

static double[][] A = new double[][] { new double[] { 1, 0, 0, 0, 0 }, new double[] { 0.760652602, 1, 0, 0, 0 }, new double[] { 0.373419404, 0.760537565, 1, 0, 0 }, new double[] { 0.136996731, 0.373331934, 0.760422587, 1, 0 }, new double[] { 0.040625222, 0.136953801, 0.373244464, 0.76030755, 1 } }; 
static double[][] b = new double[][] { new double[] { 2017159 }, new double[] { 1609660 }, new double[] { 837732.8125 }, new double[] { 330977.3125 }, new double[] { 87528.38281 } }; 

static void Main(string[] args) 
{ 
    using(Problem problem = new Problem(Solver_Type.Minimize, 5, 0)) 
    { 
     problem.VarDecision.LowerBound.Array = new double[] { 0.0, 0.0, 0.0, 0.0, 0.0 }; 
     problem.VarDecision.UpperBound.Array = new double[] { Constants.PINF, Constants.PINF, Constants.PINF, Constants.PINF, Constants.PINF }; 

     problem.Evaluators[Eval_Type.Function].OnEvaluate += new EvaluateEventHandler(SumOfSquaredErrors); 

     problem.ProblemType = Problem_Type.OptNLP; 

     problem.Solver.Optimize(); 

     Optimize_Status status = problem.Solver.OptimizeStatus; 

     Console.WriteLine(status.ToString()); 
     foreach(double x in problem.VarDecision.FinalValue.Array) 
     { 
      Console.WriteLine(x); 
     } 
    } 
} 

static Engine_Action SumOfSquaredErrors(Evaluator evaluator) 
{ 
    double[][] x = new double[evaluator.Problem.Variables[0].Value.Array.Length][]; 
    for(int i = 0; i < x.Length; i++) 
    { 
     x[i] = new double[1] { evaluator.Problem.Variables[0].Value.Array[i] }; 
    } 

    double[][] b_calculated = MatrixMultiply(A, x); 

    double sum_sq = 0.0; 
    for(int i = 0; i < b_calculated.Length; i++) 
    { 
     sum_sq += Math.Pow(b_calculated[i][0] - b[i][0], 2); 
    } 
    evaluator.Problem.FcnObjective.Value[0] = sum_sq; 

    return Engine_Action.Continue; 
} 

static double[][] MatrixMultiply(double[][] left, double[][] right) 
{ 
    if(left[0].Length != right.Length) 
    { 
     throw new ArgumentException(); 
    } 

    double[][] sum = new double[left.Length][]; 
    for(int i = sum.GetLowerBound(0); i <= sum.GetUpperBound(0); i++) 
    { 
     sum[i] = new double[right[i].Length]; 
    } 

    for(int i = 0; i < sum.Length; i++) 
    { 
     for(int j = 0; j < sum[0].Length; j++) 
     { 
      for(int k = 0; k < right.Length; k++) 
      { 
       sum[i][j] += left[i][k] * right[k][j]; 
      } 
     } 
    } 

    return sum; 
} 

no tengo ningún código para la Fundación Solver de Microsoft, porque no creo que la función objetivo puede ser escrito en una sola línea y que doesn' Permitir delegados como Solver SDK.

Answer 1

Una alternativa sería formular esto como un problema de programación lineal:

minimizar suma de los elementos en x

sujeto a Ax> = b

Esto debería ser bastante simple Formular usando Solver Foundation, basado en una de las muestras de LP.

actualización de julio 5

El enfoque anterior también se ve excesivamente complejo, pero tal vez esto se debe a la API en primera línea de Solver. Usando Fundación Solver de Microsoft, y minimizar la suma de los cuadrados de las diferencias, el siguiente programa:

private static void Main(string[] args) 
{ 
    var solver = SolverContext.GetContext(); 
    var model = solver.CreateModel(); 

    var A = new[,] 
     { 
      { 1, 0, 0, 0, 0 }, 
      { 0.760652602, 1, 0, 0, 0 }, 
      { 0.373419404, 0.760537565, 1, 0, 0 }, 
      { 0.136996731, 0.373331934, 0.760422587, 1, 0 }, 
      { 0.040625222, 0.136953801, 0.373244464, 0.76030755, 1 } 
     }; 
    var b = new[] { 2017159, 1609660, 837732.8125, 330977.3125, 87528.38281 }; 

    var n = A.GetLength(1); 
    var x = new Decision[n]; 
    for (var i = 0; i < n; ++i) 
     model.AddDecision(x[i] = new Decision(Domain.RealNonnegative, null)); 

    // START NLP SECTION 
    var m = A.GetLength(0); 
    Term goal = 0.0; 
    for (var j = 0; j < m; ++j) 
    { 
     Term Ax = 0.0; 
     for (var i = 0; i < n; ++i) Ax += A[j, i] * x[i]; 
     goal += Model.Power(Ax - b[j], 2.0); 
    } 
    model.AddGoal(null, GoalKind.Minimize, goal); 
    // END NLP SECTION 

    var solution = solver.Solve(); 
    Console.WriteLine("f = {0}", solution.Goals.First().ToDouble()); 
    for (var i = 0; i < n; ++i) Console.WriteLine("x[{0}] = {1}", i, x[i].GetDouble()); 
} 

genera la siguiente solución, que debe estar en línea con la solución de la hoja de Excel vinculada:

f = 254184688.179922 
x[0] = 2017027.31820845 
x[1] = 76226.6063397686 
x[2] = 26007.3375581303 
x[3] = 1.00650383558278E-07 
x[4] = 4.18546775823669E-09 

Si no me equivoco, a diferencia de GRG, Solver Foundation no puede admitir restricciones generales no lineales listas para usar, creo que necesitará complementos adicionales para manejarlas. Para su problema, esto no es un problema.

Para completar, a formular el problema LP lugar, intercambiar el código entre START NLP SECCIÓN y NLP FIN SECCIÓN con el siguiente código:

var m = A.GetLength(0); 
    var constraints = new Term[m]; 
    for (var j = 0; j < m; ++j) 
    { 
     Term Ax = 0.0; 
     for (var i = 0; i < n; ++i) Ax += A[j, i] * x[i]; 
     model.AddConstraint(null, constraints[j] = Model.GreaterEqual(Ax, b[j])); 
    } 
    model.AddGoal(null, GoalKind.Minimize, Model.Sum(x)); 

que producirá la siguiente salida (nota que las funciones objetivas son diferentes en los dos casos, de ahí las grandes diferencias en f):

f = 2125502.27815564 
x[0] = 2017159 
x[1] = 75302.7580022821 
x[2] = 27215.9247379241 
x[3] = 5824.5954154355 
x[4] = 0