c# - 将 MS Solver Foundation 与矩阵乘法结合使用

Question

我正在尝试使用 Microsoft Solver Foundation优化我涉及矩阵乘法的问题。我可以使用 Excel 的求解器来执行此操作，但我试图将其集成到 C# 中，但遇到了麻烦。这是一个带有示例的描述:

假设您有一个 (3x3) 矩阵 y，定义为:

Double[][] y = 
{
    new Double[] { 5, 1, 0 },
    new Double[] { 1, 9, 1 },
    new Double[] { 0, 1, 9 },
};

我想找到 (1x3) 矩阵 x 使得:x * y * x' 最小化。 此外，x 值的总和必须为 1，并且所有 x 值都不能小于 0。

这是我目前的代码:

 SolverContext context = SolverContext.GetContext();             // Get context environment
 Model model = context.CreateModel();                            // Create a new model

 Decision d1 = new Decision(Domain.RealNonnegative, "d1");       // First item in "x" vector (must be >= 0)
 Decision d2 = new Decision(Domain.RealNonnegative, "d2");       // Second item in "x" vector (must be >= 0)
 Decision d3 = new Decision(Domain.RealNonnegative, "d3");       // Third item in "x" vector (must be >= 0)
 model.AddDecisions(d1, d2, d3);                                 // Add these to the model (this is where the outputs will be stored)

 model.AddConstraints("limits",                                  // Add constraints
     0 <= d1 <= 1,                                               // Each item must be between 0 and 1
     0 <= d2 <= 1,
     0 <= d3 <= 1,
     d1 + d2 + d3 == 1);                                         // All items must add up to 1

我坚持的部分是当你告诉它你想要最小化什么时:

 model.AddGoal("min", GoalKind.Minimize, /* What goes here? */);

这部分通常包含一个方程(例如 d1 * d2 + d3)，但矩阵乘法并不那么简单。

我可以创建一个执行乘法并返回 double 的函数，但是 AddGoal() 需要一个 Term对象，而且我还必须对 Decision 对象进行算术运算。

或者，我可以将这个乘法分解为一个巨大的 string 表达式(我已经这样做了)，但如果我不必这样做，我会更愿意。 (这个字符串看起来像:"d1 * 5 + d2 * 1 + d3 * 0 ...")

有什么想法吗？ 谢谢。

PS:正确答案(根据 Excel)是:

d1 = 0.503497
d2 = 0.216783
d3 = 0.27972

注意:解决方案必须可扩展以具有 n 个“决策”

Answer 1

以下产生预期的解决方案:

using System;
using Microsoft.SolverFoundation.Services;

namespace akMSFStackOverflow
{
    class Program
    {
        static void Main(string[] args)
        {
            Double[,] y = 
            {
                { 5, 1, 0 },
                { 1, 9, 1 },
                { 0, 1, 9 },
            };

            Term goal;
            Term[,] tx;
            Term[,] ty;

            SolverContext context = SolverContext.GetContext();             // Get context environment
            Model model = context.CreateModel();                            // Create a new model

            Decision d1 = new Decision(Domain.RealNonnegative, "d1");       // First item in "x" vector (must be >= 0)
            Decision d2 = new Decision(Domain.RealNonnegative, "d2");       // Second item in "x" vector (must be >= 0)
            Decision d3 = new Decision(Domain.RealNonnegative, "d3");       // Third item in "x" vector (must be >= 0)
            model.AddDecisions(d1, d2, d3);                                 // Add these to the model (this is where the outputs will be stored)

            model.AddConstraints("limits",                                  // Add constraints
                0 <= d1 <= 1,                                               // Each item must be between 0 and 1
                0 <= d2 <= 1,
                0 <= d3 <= 1,
                d1 + d2 + d3 == 1);                                         // All items must add up to 1

            ty = matrix(y);
            tx = new Term[,] { { d1, d2, d3 } };

            goal = matMult(matMult(tx, ty), transpose(tx))[0, 0];

            model.AddGoal("goal", GoalKind.Minimize, goal);

            // Specifying the IPM solver, as we have a quadratic goal 
            Solution solution = context.Solve(new InteriorPointMethodDirective());


            Report report = solution.GetReport();
            Console.WriteLine("x {{{0}, {1}, {2}}}", d1, d2, d3);
            Console.Write("{0}", report); 

        }


        static Term[,] matrix(Double[,] m)
        {
            int rows = m.GetLength(0);
            int cols = m.GetLength(1);
            Term[,] r = new Term[rows, cols];

            for (int row = 0; row < rows; row++)
                for (int col = 0; col < cols; col++)
                    r[row, col] = m[row, col];

            return r;
        }

        static Term[,] matMult(Term[,] a, Term[,] b)
        {
            int rows = a.GetLength(0);
            int cols = b.GetLength(1);
            Term[,] r = new Term[rows, cols];

            for (int row = 0; row < rows; row++)
                for (int col = 0; col < cols; col++)
                {
                    r[row,col] = 0;
                    for (int k = 0; k < a.GetLength(1); k++)
                    {
                        r[row, col] += a[row, k] * b[k, col];
                    }
                }

            return r;
        }

        static Term[,] transpose(Term[,] m)
        {
            int rows = m.GetLength(0);
            int cols = m.GetLength(1);
            Term[,] r = new Term[cols, rows];

            for (int row = 0; row < rows; row++)
                for (int col = 0; col < cols; col++)
                {
                    r[col, row] = m[row, col];
                }

            return r;
        }
    }
}