Nathan Brixius

Erwin, a modeling consultant and top Solver Foundation user, encountered some problems trying to do two-way data binding using DataTable objects.  There are more details on this discussion thread.  Ross, a member of the Solver Foundation team, was kind enough to code up a workaround for Erwin's example.  In addition to this CS file, you will need to create a new DBML file called SampleDataContext, as I described in a previous post.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Data;
using System.Data.OleDb;
using System.Data.Linq;
using System.Text;
using Microsoft.SolverFoundation.Services;
using System.IO;
 
namespace OML1
{
    class Test
    {
        static void Main(string[] args)
        {
            Test t = new Test();
            t.Solve();
        }
 
        // Holds the OML model
        string strModel = @"Model[
              Parameters[Sets,I],
              Parameters[Reals,p[I]],
 
              Decisions[Reals[0,Infinity],x[I]],
 
              Constraints[
                 Foreach[{i,I}, x[i]==p[i]]
              ]
           ]";
 
 
        //  SFS
        SolverContext context;
        SampleDataContext data;
 
 
        //  Constructor
        public Test()
        {
            context = SolverContext.GetContext();
            data = new SampleDataContext("Data Source=Sql_server_name;Initial Catalog=DataPartitionAllocation20_5;Integrated Security=True");
            context.DataSource = data;
        }
 
 
        // Solve the problem
        public void Solve()
        {
            context.LoadModel(FileFormat.OML, new StringReader(strModel));
 
            Parameter p = context.CurrentModel.Parameters.First(q => q.Name == "p");
            p.SetBinding(data.P, "value", new string[] { "index" });
 
            Decision x = context.CurrentModel.Decisions.First(d => d.Name == "x");
            x.SetBinding(data.X, "value", new string[] { "index" }); 
 
            Solution solution = context.Solve();
            Console.Write("{0}", solution.GetReport());
 
            context.PropagateDecisions();
 
        }
 
    }
 
}

Here is my as-promised "part 2" to my traveling salesman sample. My goal is not to try and do any "production level" modeling, but instead "point the way" by highlighting different aspects of Solver Foundation's architecture. This time I want to talk about directives, plug-in solvers, and model reusability. The Solver Foundation architecture provides the ability to plug-in other solvers to handle one or more problem types. As of version 1.1, the Gurobi solver is Solver Foundation's default MIP (mixed integer programming) solver. Gurobi is a high-performance solver which is capable of taking on industrial strength problems, so my piddling little TSP example is no match for it. Gurobi has a number of parameters that can be set to tune solver performance - you set them by means of Directives. First, let's modify our code to turn on logging. Change the context.Solve() call to:

      GurobiDirective gurobiDirective = new GurobiDirective();
      gurobiDirective.OutputFlag = true;
      Solution solution = context.Solve(gurobiDirective);

Each plug-in solver also defines its own directives - therefore the full power and flexibility of plug-in solvers are available to you as a programmer. In this case I have just flipped on the logging bit - here's what I see.

Optimize a model with 185 Rows, 380 Columns and 1199 NonZeroes
Presolve removed 1 rows and 185 columns
Presolve time: 0.30s
Presolved: 184 Rows, 195 Columns, 832 Nonzeros
Objective GCD is 1

Root relaxation: objective 2.783286e+03, 36 iterations, 0.02 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0  2783.2857    0   14          -  2783.2857     -      -    0s
     0     0  3006.1964    0   29          -  3006.1964     -      -    0s
     0     0  3006.1964    0   29          -  3006.1964     -      -    0s
     0     0  3006.1964    0   25          -  3006.1964     -      -    0s
     0     0  3006.1964    0   25          -  3006.1964     -      -    0s
H    0     0                       6029.0000  3006.1964  50.1%     -    0s
     0     0  3006.1964    0   25  6029.0000  3006.1964  50.1%     -    0s
H    0     2                       3805.0000  3006.1964  21.0%     -    0s
     0     2  3006.1964    0   25  3805.0000  3006.1964  21.0%     -    0s
*   75    54              29       3644.0000  3077.0000  15.6%   8.2    0s
*  129    80              27       3381.0000  3077.0000  8.99%   8.7    0s
*  158    62              30       3323.0000  3077.0000  7.40%   8.1    0s

Cutting planes:
  Learned: 3
  Gomory: 1
  MIR: 3

Explored 1231 nodes (5693 simplex iterations) in 0.69 seconds
Thread count was 2 (of 2 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective 3.3230000000e+03, best bound 3.3230000000e+03, gap 0.0000%

For large models it is sometimes useful to tweak the solver options for performance. For example, if I wanted to spend a slightly larger amount of time doing MIP heuristics, and I wanted to eliminate presolve, then I can set those options easily:

      gurobiDirective.Heuristics = 0.1;
      gurobiDirective.Presolve = PresolveLevel.None;

There is a full reference in the Gurobi Solver Programming Primer that ships with Solver Foundation. My point is that you do not get a "least common denominator" solution - you get all the bells and whistles.

Gurobi is just one example - there are a wide range of plug-ins available, covering most of the popular LP and MIP solvers - check out solverfoundation.com for the full list. To use a plug-in solver, you need to: 1) install the solver (duh), 2) change your exe.config or web.config to point to the solver. No code changes are required. You can find additional documentation on the default MIP directives and the third-party Solver configuration in the SFS Programming Primer.

Last thing - modern programming languages like C# are cool because they promote reusability and maintainability. Don't forget about that when you are writing SFS code. It's common for models in a particular vertical (say, finance or transportation) to have common "blocks" of goals or constraints that are re-used over and over again. It's pretty easy to build reusable model libraries using extension methods. For example, take the assignment constraints in my TSP example. I could factor them out:

  public static class ModelingExtensions {
    public static void AssignmentConstraintsNoDiag(this Model model, Set s, Decision assign) {
      model.AddConstraint("A1", Model.ForEach(s, i => Model.Sum(Model.ForEachWhere(s, j => assign[i, j], j => i != j)) == 1));
      model.AddConstraint("A2", Model.ForEach(s, j => Model.Sum(Model.ForEachWhere(s, i => assign[i, j], i => i != j)) == 1));
    }
  }

Now I can just say model.AssignmentConstraintsNoDiag(city, assign). I could add other variations inside of ModelingExtensions as well. I think this aspect of Solver Foundation has not gotten enough attention - maybe because it doesn't make for a cool demo!

On the Solver Foundation MSDN forum there was a question about how to read model data from a DB and use it within a Solver Foundation model.  In this post I will extend my production planning sample to use LINQ to SQL.  To follow along at home you will need to have a recent version of SQL Server installed locally, and some basic knowledge of how to create SQL tables.  You should also have compiled and run the code from my previous post.

Step 1: Create and populate the DB 

The first step is to create tables corresponding to the entities in my model.  I created a very simple DB with three tables: Countries, Products, and Yields.  The Yields table has foreign key constraints to the Countries and Products tables.  Here is a diagram:

To populate the DB I just wrote a script that inserts my problem data, and ran itin SQL Management Studio.  Here's the script (and forgive my SQL):

GO 

DELETE FROM Yields
DELETE FROM Products
DELETE FROM Countries
GO

INSERT INTO Countries (Id, Name, Limit, Cost)
VALUES (0, 'SA', 9000, 20)

INSERT INTO Countries (Id, Name, Limit, Cost)
VALUES (1, 'VZ', 6000, 15)

GO

INSERT INTO Products (Id, Name, Demand)
VALUES (0, 'Gas', 1900)

INSERT INTO Products (Id, Name, Demand)
VALUES (1, 'Jet Fuel', 1500)

INSERT INTO Products (Id, Name, Demand)
VALUES (2, 'Lubricant', 500)

GO

INSERT INTO Yields (CountryId, ProductId, Value)
VALUES (0, 0, 0.3)
INSERT INTO Yields (CountryId, ProductId, Value)
VALUES (1, 0, 0.4)

INSERT INTO Yields (CountryId, ProductId, Value)
VALUES (0, 1, 0.4)
INSERT INTO Yields (CountryId, ProductId, Value)
VALUES (1, 1, 0.2)

INSERT INTO Yields (CountryId, ProductId, Value)
VALUES (0, 2, 0.2)
INSERT INTO Yields (CountryId, ProductId, Value)
VALUES (1, 2, 0.3)

 

Step 2: Create Entity and DataContext classes in Visual Studio

Scott Guthrie's blog (and the MSDN docs) show you exactly how to do this: 

1. Add a new "Linq to SQL file" to your project called Petrochem.dbml.
2. Bring up the Server Explorer window, connect to your database and drag the tables into the dbml window.  Visual Studio will automatically create a Datacontext class (mine is called PetrochemDataContext) and an entity class for each table that you include.

Step 3: Modify Solver Foundation Services data binding code

This is in fact very easy because Solver Foundation Services was designed to work well with LINQ.  Take the PetrochemDataBinding sample from last time, and change the SetBinding statements to work with the PetrochemDataContext class instead of a hardcoded DataSet.  The code is almost identical:

  private static void PetrochemLinqDataBinding() {
      SolverContext context = SolverContext.GetContext();
      context.ClearModel();
      Model model = context.CreateModel();

      PetrochemDataContext db = new PetrochemDataContext();

      Set products = new Set(Domain.Any, "products");
      Set countries = new Set(Domain.Any, "countries");

      Parameter demand = new Parameter(Domain.Real, "demand", products);
      demand.SetBinding(db.Products, "Demand", "Id");

      Parameter yield = new Parameter(Domain.Real, "yield", products, countries);
      yield.SetBinding(db.Yields, "Value", "ProductId", "CountryId");

      Parameter limit = new Parameter(Domain.Real, "limit", countries);
      limit.SetBinding(db.Countries, "Limit", "Id");

      Parameter cost = new Parameter(Domain.Real, "cost", countries);
      cost.SetBinding(db.Countries, "Cost", "Id");

      model.AddParameters(demand, yield, limit, cost);

      Decision produce = new Decision(Domain.RealNonnegative, "produce", countries);
      model.AddDecision(produce);

      model.AddGoal("goal", GoalKind.Minimize, Model.Sum(Model.ForEach(countries, c => cost[c] * produce[c])));

      model.AddConstraint("Demand",
        Model.ForEach(products, p => Model.Sum(Model.ForEach(countries, c => yield[p, c] * produce[c])) >= demand[p])
        );

      model.AddConstraint("Production limit",
        Model.ForEach(countries, c => produce[c] <= limit[c])
        );

      Solution solution = context.Solve(new SimplexDirective());
      Report report = solution.GetReport();
      Console.WriteLine(report);
    }

That's all there is to it!  Note that instead of passing the entire collection (e.g. db.Countries) you could easily use LINQ statements or stored procedures, or whatever you like.

Here's an example that I walked through during yesterday's INFORMS session.  Erwin has two blog postings about Solver Foundation and the traveling salesman problem, but I want to throw in my two cents because I want to emphasize a couple of points:

1. By combining C# and Solver Foundation Services it is possible to express complex models clearly and succinctly.
2. It is very easy to build powerful, reusable model libraries using C# and Solver Foundation Services.
3. Solver Foundation Services code can be used in many different application environments (ASP.Net, silverlight, DB, command line apps, WPF, …) with minimal changes.

The traveling salesman problem is a classical problem in computer science, and you should bow your head in shame if you don't know about it (and turn in your conference badge if you happen to be in Phoenix).  A salesperson needs to make a tour of a number of cities.  The restrictions are that she wants to visit each city once and only once, and she wants to minimize the distance travelled.  This is perhaps the definitive example of an NP-hard problem.

TSP can be solved using mixed integer programming - optimizing a linear goal with linear constraints, where some of the decision variables are integer.  In this first post I will show how to formulate and solve a TSP model using Solver Foundation Services.  In my second post I will show how to use the Gurobi MIP solver using SFS.   There are many different ways to model the TSP - here is a nice introduction.  My goal is to provide a clear, complete example - not build a "production level" TSP model, so I am going to choose a model formulation that dates back to 1960!  First, I need to establish a couple of building blocks that will help me construct the data for the model.  We need to know the distances between each pair of cities.  Typically we are provided the coordinates of the cities and need to derive the distances.  So I will introduce a Coordinate class that contains properties for the (x, y) coordinates, and properties to convert to latitude and longitude.  Finally, a method that computes the distance between points.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Microsoft.SolverFoundation.Services;
using SolverFoundation.Plugin.Gurobi;

namespace Microsoft.SolverFoundation.Samples.TravelingSalesman {
  class Program {
    // TSP coordinate.
    public class Coordinate {
      public int Name { get; set; }

      // X-coordinate (from TSPLIB)
      public double X { get; set; }

      // Y-coordinate (from TSPLIB)
      public double Y { get; set; }

      public Coordinate(int name, double x, double y) {
        Name = name;
        X = x;
        Y = y;
      }

      // Latitude in radians.
      public double Latitude {
        get { return Math.PI * (Math.Truncate(X) + 5 * (X - Math.Truncate(X)) / 3) / 180; }
      }

      // Longitude in radians.
      public double Longitude {
        get { return Math.PI * (Math.Truncate(Y) + 5 * (Y - Math.Truncate(Y)) / 3) / 180; }
      }

      // Geographic distance between two points (as an integer).
      public int Distance(Coordinate p) {
        double q1 = Math.Cos(Longitude - p.Longitude);
        double q2 = Math.Cos(Latitude - p.Latitude);
        double q3 = Math.Cos(Latitude + p.Latitude);
        // There may rounding difficulties her if the points are close together...just sayin'.
        return (int)(6378.388 * Math.Acos(0.5 * ((1 + q1) * q2 - (1 - q1) * q3)) + 1);
      }
    }

    // TSP city-city arc.
    // 
    public class Arc {
      public int City1 { get; set; }
      public int City2 { get; set; }
      public double Distance { get; set; }
    }

    // Burma14 from TSPLIB. Optimal tour = 3323.
    private static Coordinate[] data = new Coordinate[] {
      new Coordinate(0, 16.47, 96.10),
      new Coordinate(1, 16.47, 94.44),
      new Coordinate(2, 20.09, 92.54),
      new Coordinate(3, 22.39, 93.37),
      new Coordinate(4, 25.23, 97.24),
      new Coordinate(5, 22.00, 96.05),
      new Coordinate(6, 20.47, 97.02),
      new Coordinate(7, 17.20, 96.29),
      new Coordinate(8, 16.30, 97.38),
      new Coordinate(9, 14.05, 98.12),
      new Coordinate(10, 16.53, 97.38),
      new Coordinate(11, 21.52, 95.59),
      new Coordinate(12, 19.41, 97.13),
      new Coordinate(13, 20.09, 94.55)
    };

(The data for this 14-city problem comes from the TSPLIB library). If you've been following my blog you know that the building blocks of a Solver Foundation model are: sets, parameters, decisions, goals, and constraints. I am going to implement a simple formulation that is centered around the following (indexed) decisions:

• Assign[i,j]: this is equal to 1 if the optimal tour contains a trip (or arc) from city i to city j.
• Rank[i]: this is equal to the number of cities visited after arriving at city i. 

We have one parameter in our model:

• Distance[I,j]: the distance from city i to city j. 

With that in mind, here's the model.  Explanation of the goals and constraints follow.

    static void Main(string[] args) {
      SolverContext context = SolverContext.GetContext();
      Model model = context.CreateModel();

      // ------------
      // Parameters
      Set city = new Set(Domain.IntegerNonnegative, "city");
      Parameter dist = new Parameter(Domain.Real, "dist", city, city);
      var arcs = from p1 in data
                 from p2 in data
                 select new Arc { City1 = p1.Name, City2 = p2.Name, Distance = p1.Distance(p2) };
      dist.SetBinding(arcs, "Distance", "City1", "City2");
      model.AddParameters(dist);

      // ------------
      // Decisions
      Decision assign = new Decision(Domain.IntegerRange(0, 1), "assign", city, city);
      Decision rank = new Decision(Domain.RealNonnegative, "rank", city);
      model.AddDecisions(assign, rank);

      // ------------
      // Goal: minimize the length of the tour.
      Goal goal = model.AddGoal("TourLength", GoalKind.Minimize,
        Model.Sum(Model.ForEach(city, i => Model.ForEachWhere(city, j => dist[i, j] * assign[i, j], j => i != j))));

      // ------------
      // Enter and leave each city only once.
      int N = data.Length;
      model.AddConstraint("assign 1",
        Model.ForEach(city, i => Model.Sum(Model.ForEachWhere(city, j => assign[i, j],
          j => i != j)) == 1));
      model.AddConstraint("assign 2",
        Model.ForEach(city, j => Model.Sum(Model.ForEachWhere(city, i => assign[i, j], i => i != j)) == 1));
      model.AssignmentNoDiag(city, assign);

      // Forbid subtours (Miller, Tucker, Zemlin - 1960...)
      model.AddConstraint("no subtours",
        Model.ForEach(city,
          i => Model.ForEachWhere(city,
            j => rank[i] + 1 <= rank[j] + N * (1 - assign[i, j]),
            j => Model.And(i != j, i >= 1, j >= 1)
          )
        )
      );

      Solution solution = context.Solve();

      // Retrieve solution information.
      Console.WriteLine("Cost = {0}", goal.ToDouble());
      Console.WriteLine("Tour:");
      var tour = from p in assign.GetValues() where (double)p[0] > 0.9 select p[2];
      foreach (var i in tour.ToArray()) {
        Console.Write(i + " -> ");
      }
      Console.WriteLine();
    }
  }

In my humble opinion, the "Parameter data =" line is an awesome example of the power of LINQ data binding in Solver Foundation.  We generate the 2D matrix of distances using a single LINQ expression.  It would be incredibly easy to change the code to retrieve the coordinate data from a database (perhaps using a LINQ expression once again), a file, or even a user application.

The goal is straightforward: minimize the distance traveled.  This is a product of the selected arcs and the distance matrix.   We have two types of constraints:

• Assignment constraints: these ensure that we enter and leave each city only once.

• Subtour constraints: these ensure that we do not have any subtours.  In a four city problem {A, B, C, D}, for example, we cannot have two cycles (A, B), (C, D).  We need to have one tour that contains all the cities.

The assignment constraints are easy using the ForEach and ForEachWhere operations.  I use ForEachWhere because I want to disallow arcs that enter and leave the same city - that doesn't make sense.  The subtour constraint is a little more complicated.   It relates the "assign" and "rank" decisions.  The key fact is that if there is an arc from city i to city j, rank[i] + 1 == j.  Of course, if the (i, j) arc is not part of the optimal tour then all bets are off.   Last note: notice that I can mix parameters, decisions, and C# variables in my expressions. 

 
Once we solve the model we want to print out the cost, and the optimal tour.  Getting the cost is very easy using goal.ToDouble().  We can get the tour using either Assign or Rank.  I have chosen to use Assign because it gives me another opportunity to use LINQ.  When you call GetValues() on a decision, you get arrays that contain the value along with the indexes for each decision.  In this case, the last entry in the array is the one we are interested in.  There are other ways to conveniently query decsision results, I'll save that for another time.
 
The next post will show how we can use Solver Foundation's plug-in model to tune the behavior of the Gurobi MIP solver.
 

In this post I am going to present two complete C# programs for modeling and solving a simple production problem using Solver Foundation. In this example we have two refineries (located in Saudi Arabia and Venzuela) that produce three products: gasoline, jet fuel, and lubricant. The goal is to minimize production costs, which depend on location. There is demand for each product which must be met. Finally, each production site has a limited production capacity.

The code for the simple case is below.  Creating a model amounts to adding the goals, decisions, and constraints using the appropriate Add method.  The signature for each is similar: the first argument is the name and the second argument is a term.  Terms are created by combining decisions or parameters using the usual operators, or by using static methods on the Model class (as we will see in our second example).  Finally, notice the call to SolverContext.Solve().  I have supplied a Directive that specifies that the Simplex solver should be used.

using System;
using System.Collections.Generic;
using System.Data;
using System.Linq;
using System.Text;
using Microsoft.SolverFoundation.Services;

namespace Microsoft.SolverFoundation.Samples.Petrochem {
  class Program {
    static void Main(string[] args) {
      PetrochemSimple();
    }

    private static void PetrochemSimple() {
      SolverContext context = SolverContext.GetContext();
      context.ClearModel();
      Model model = context.CreateModel();

      Decision sa = new Decision(Domain.RealRange(0, 9000), "SA");
      Decision vz = new Decision(Domain.RealRange(0, 6000), "VZ");
      model.AddDecisions(sa, vz);

      model.AddGoal("goal", GoalKind.Minimize, 20 * sa + 15 * vz);

      model.AddConstraint("demand1", 0.3 * sa + 0.4 * vz >= 1900);
      model.AddConstraint("demand2", 0.4 * sa + 0.2 * vz >= 1500);
      model.AddConstraint("demand3", 0.2 * sa + 0.3 * vz >= 500);

      Solution solution = context.Solve(new SimplexDirective());
      Report report = solution.GetReport();
      Console.WriteLine(report);
    }
  }

Now let's reimplement the same example using data binding. Data binding is a powerful mechanism for creating large, maintainable models. Notice that in my first example, the numeric data such as the yields, the demands, and capacities were expressed directly in the terms. The first step in using Solver Foundation data binding is to lift these values into Parameters. It is often useful to create indexed parameters using Sets. In this example, there are two clearly defined Sets: the set of countries, and the set of products. In my example below I create a DataSet which contains the data for each of my parameters. (The GetData() method is just an example, it's not pretty but it is needed to complete the example.) Then I create a series of indexed parameters. For each of them I call the SetBinding method to associate the data with the parameter. In addition to the data, SetBinding also requires the caller to indicate which property specifies the values for the parameter. If the parameter is indexed, I also need to specify the parameters of the index properties. Since I am working with DataTables, these are simply column names. Notice that I could swap in any other data source that is enumerable - in particular LINQ works really well with Solver Foundation parameters.

After the parameters are created, I define the decisions, goals, and constraints. Notice that there is only one decision - it is indexed. The Model.Sum and Model.Foreach operations allow me to define a series of constraints over one or more indexed sets in one single statement. This means that if I were to add more countries or products, my model definition would not change at all.

    private static void PetrochemDataBinding() {
      SolverContext context = SolverContext.GetContext();
      context.ClearModel(); 
      Model model = context.CreateModel();

      // Retrieve the problem data.
      DataSet data = GetData();

      Set products = new Set(Domain.Any, "products");
      Set countries = new Set(Domain.Any, "countries");
      
      Parameter demand = new Parameter(Domain.Real, "demand", products);
      demand.SetBinding(data.Tables["Demand"].AsEnumerable(), "Demand", "Product");

      Parameter yield = new Parameter(Domain.Real, "yield", products, countries);
      yield.SetBinding(data.Tables["Yield"].AsEnumerable(), "Yield", "Product", "Country");

      Parameter limit = new Parameter(Domain.Real, "limit", countries);
      limit.SetBinding(data.Tables["Limit"].AsEnumerable(), "Limit", "Country");

      Parameter cost = new Parameter(Domain.Real, "cost", countries);
      cost.SetBinding(data.Tables["Cost"].AsEnumerable(), "Cost", "Country");

      model.AddParameters(demand, yield, limit, cost);

      Decision produce = new Decision(Domain.RealNonnegative, "produce", countries);
      model.AddDecision(produce);

      model.AddGoal("goal", GoalKind.Minimize, Model.Sum(Model.ForEach(countries, c => cost[c] * produce[c])));

      model.AddConstraint("Demand", 
        Model.ForEach(products, p => Model.Sum(Model.ForEach(countries, c => yield[p, c] * produce[c])) >= demand[p])
        );

      model.AddConstraint("Production limit",
        Model.ForEach(countries, c => produce[c] <= limit[c])
        );

      Solution solution = context.Solve(new SimplexDirective());
      Report report = solution.GetReport();
      Console.WriteLine(report);
    }

    private static DataSet GetData() {
      string[] products = new string[] { "Gas", "Jet Fuel", "Lubricant" };
      string[] countries = new string[] { "SA", "VZ" };

      double[][] yield = new double[][] {
        new double[] { 0.3, 0.4 }, 
        new double[] { 0.4, 0.2 }, 
        new double[] { 0.2, 0.3 } 
      };
      double[] demand = new double[] { 1900, 1500, 500 };
      double[] limit = new double[] { 9000, 6000 };
      double[] cost = new double[] { 20, 15 };

      DataSet dataSet = new DataSet();
      #region Fill DataSet
      DataTable table = new DataTable("Yield");
      dataSet.Tables.Add(table);
      table.Columns.Add("Product", typeof(string));
      table.Columns.Add("Country", typeof(string));
      table.Columns.Add("Yield", typeof(double));

      for (int p = 0; p < products.Length; p++) {
        for (int c = 0; c < countries.Length; c++) {
          DataRow row = table.NewRow();
          row[0] = products[p];
          row[1] = countries[c];
          row[2] = yield[p][c];
          table.Rows.Add(row);
        }
      }

      table = new DataTable("Demand");
      dataSet.Tables.Add(table);
      table.Columns.Add("Product", typeof(string));
      table.Columns.Add("Demand", typeof(double));
      for (int p = 0; p < products.Length; p++) {
        DataRow row = table.NewRow();
        row[0] = products[p];
        row[1] = demand[p];
        table.Rows.Add(row);
      }

      table = new DataTable("Limit");
      dataSet.Tables.Add(table);
      table.Columns.Add("Country", typeof(string));
      table.Columns.Add("Limit", typeof(double));
      for (int c = 0; c < countries.Length; c++) {
        DataRow row = table.NewRow();
        row[0] = countries[c];
        row[1] = limit[c];
        table.Rows.Add(row);
      }

      table = new DataTable("Cost");
      dataSet.Tables.Add(table);
      table.Columns.Add("Country", typeof(string));
      table.Columns.Add("Cost", typeof(double));
      for (int c = 0; c < countries.Length; c++) {
        DataRow row = table.NewRow();
        row[0] = countries[c];
        row[1] = cost[c];
        table.Rows.Add(row);
      }
      #endregion

      return dataSet;
    }

A few of us on the Solver Foundation team will be presenting at the INFORMS Practice Conference in Phoenix, Arizona.  I'm looking forward to it!  We have both a workshop and a session - I will also be participating in Gurobi's session to show off its integration with Solver Foundation.  Here are the details - looking forward to seeing you there.

Technology Workshop
Sunday, April 26, 9AM - 12PM
Microsoft Solver Foundation is a pure, managed code runtime for mathematical programming, modeling, and optimization. Microsoft Solver Foundation provides solvers and services to a broad community of users: from Excel users and analysts to programmers working on business critical scheduling, configuration, risk management and planning solutions. It provides services for model validation, parallel solving and workload scheduling, model interchange, and declarative data binding via LINQ and other NETFx technologies. As an open framework designed for third party extensibility, it exposes facilities for users to plug-in their own solvers while still leveraging all of the modeling services and capabilities of Microsoft Solver Foundation.

Software Tutorials
Monday, April 27, 9:10AM - 10:00AM
Solver Foundation is a framework and managed code runtime for mathematical programming, modeling and optimization. This tutorial will focus on the technical overview of Solver Foundation, and how Third-party solver vendors, modelers and solution specialists can leverage Solver Foundation from all CLS-compliant languages and Microsoft Office.

I am pleased to announce that Solver Foundation v1.2 is live on solverfoundation.com! Our goal in 1.2 was to make a few key improvements to address feedback that we have received from partners and customers. If you have feedback, questions, suggestions…please post it here, or on our MSDN forum. Positive or negative - that's okay. I am particularly interested in hearing how you use Solver Foundation in your application - or what features you would like to see.

Back to 1.2:

• The popular FICO Xpress is now a "certified partner". That means if you are an Xpress user, you can use Solver Foundation's plug-in model to write .Net applications without a hitch.
• Amsterdam Optimization is a certified modeling consultant for Solver Foundation. I've linked to Erwin's blog before - he's a member of this group. We are happy to be partnering with them.
• MIP start functionality has been added. This allows you to give a starting solution to a MIP model using SetInitialValue on Decision.
• Bug fixes - in particular we've taken care of a few QP issues.

I have been brushing up on keyboard shortcuts ever since I installed the Windows 7 beta on my laptop.  Tim Sneath has a great summary of keyboard shortcuts here, in particular I LOVE the window management shortcuts.  If you have IE8 you also probably know that Ctrl+Tab switches between browser tabs.  What I did not know, and just discovered by accident, is that Ctrl+Tab works in Visual Studio too.  You can use it to switch between active files and tool windows (e.g. the Solution Explorer or Class View).  The arrow keys work too.  I tend to have lots of active files at once, so this is a nice timesaver.

Feel free to share other cool shortcuts in the comments.

A side note: I am a Unix guy by background, and during my first two years at Microsoft I did all my development using the command shell + emacs.  But given the improvements that Visual Studio has made over the last 5 years, there is no way I would ever go back now!  I am anxiously awaiting their next release.

Two weeks back I posted two articles showing how easy it is to model critical path scheduling using Microsoft Solver Foundation.  I received a few emails asking about various extensions; I will be covering those in upcoming posts.  Julian just wrote a great blog post that covers the most commonly requested extension - resource constrained scheduling.  If you are getting started with Solver Foundation and want to see an interesting, instructive example, I encourage you to check out his post.  Two things that I really like about his OML:

1. He separates the data from the model description using Parameters.
2. He relies on the Foreach construct to define his constraints.  It results in a very clean model definition.

Julian mentions it in his post, but I also want to call out another great resource for learning OML.  Erwin Kalvalagen has written an OML tutorial which includes several interesting examples.  It's a great complement to the Excel Programming Primer that is part of the Solver Foundation documentation.

I am going to be rolling out the rest of my branch-and-bound algorithm in the next few posts.  To make that easier, in this post I introduce some common matrix, vector, and permutation methods.  It turns out that for technical computing applications, C#'s extension methods (introduced in 3.0) are awesome.  With vectors it's great because you retain the control of having direct array access, but you get the nice object-oriented notation.

I've left out comments and error-checking, and it's not super-optimized, but you get the idea.  None of these methods end up being the bottleneck.  Most of the methods are obvious, but I do want to comment on the permutation methods.  The goal of QAP is to find an optimal assignment of facilities to locations, represented as a permutation.  In the course of our branch-and-bound algorithm, we'll be making partial assignments - that is, only some of the entries in the permutation will be filled in.  By convention, p[i] == -1 will mean that facility i is unassigned.  An index where p[i] < 0 is called "unused".  When we branch, we'll want to pick an unused facility (or location), and try assigning all unused locations (or facilities) to it.   The use of C#'s iterator and yield concepts really help here.

That said, here's the code.

  public static class MatrixUtilities {
    #region Vector
    public static void ConstantFill(this T[] data, T val) {
      for (int i = 0; i < data.Length; i++) {
        data[i] = val;
      }
    }

    public static int ArgMin(this double[] data, out double best) {
      if (data.Length == 0) {
        best = Double.MinValue;
        return -1;
      }
      int iBest = 0;
      best = data[0];
      for (int i = 1; i < data.Length; i++) {
        if (best > data[i]) {
          best = data[i];
          iBest = i;
        }
      }
      return iBest;
    }

    public static int ArgMax(this double[] data, out double best) {
      if (data.Length == 0) {
        best = Double.MinValue;
        return -1;
      }
      int iBest = 0;
      best = data[0];
      for (int i = 1; i < data.Length; i++) {
        if (best < data[i]) {
          best = data[i];
          iBest = i;
        }
      }
      return iBest;
    }
    #endregion

    #region Permutation
    public static void Swap(this int[] p, int i, int j) {
      int temp = p[i];
      p[i] = p[j];
      p[j] = temp;
    }

    public static int FindUnused(this int[] data, int index) {
      foreach (int unused in data.UnusedIndices()) {
        if (index-- <= 0) {
          return unused;
        }
      }
      return -1;
    }

    public static IEnumerable UnusedIndices(this int[] data) {
      for (int i = 0; i < data.Length; i++) {
        if (data[i] < 0) {
          yield return i;
        }
      }
    }

    public static string PermutationToString(this int[] p, bool oneBased) {
      StringBuilder build = new StringBuilder(p.Length * 4);
      int width = ((int)Math.Log10(p.Length)) + 2;
      build.Append("[");
      for (int i = 0; i < p.Length; i++) {
        int index = oneBased ? p[i] + 1 : p[i];
        build.Append(index.ToString().PadLeft(width));
      }
      build.Append("]");
      return build.ToString();
    }
    #endregion

    #region Matrix

    public static string MatrixToString(this double[][] A) {
      if (A != null) {
        StringBuilder build = new StringBuilder(A.Length * A.Length * 3);
        for (int i = 0; i < A.Length; i++) {
          if (A[i] != null) {
            for (int j = 0; j < A[i].Length; j++) {
              build.AppendFormat("{0,4}", A[i][j]);
              build.Append(" ");
            }
            build.AppendLine();
          }
        }
        return build.ToString();
      }
      return null;
    }

    public static double[][] NewMatrix(int m, int n) {
      double[][] M = new double[m][];
      for (int i = 0; i < M.Length; i++) {
        M[i] = new double[n];
      }
      return M;
    }

    public static double[][] NewMatrix(int n) {
      return NewMatrix(n, n);
    }
    #endregion
  }

This is part 4 in a series of posts laying out a simple branch-and-bound solver for QAP in C#. Last time (several months ago!) I provided a simple bounding procedure for QAP.  I want to take a step back and give the high-level algorithm, then in subsequent posts I will fill in the implementation.

Branch-and-bound algorithms are a general framework for solving combinatorial optimization problems.  The idea is to repeatedly subdivide the search space into smaller subproblems of the same time.  For each subproblem we compute bounds, which may allow us to determine that the subproblem cannot lead to an improved solution to the original problem.

  /// 
Solves QAP using branch-and-bound.
  /// 
  public class BranchAndBound {

    private readonly Qap _qap;
    private BranchAndBoundResults _results;

    /// 
Create a new instance.
    /// 
    public BranchAndBound(Qap qap) {
      _qap = qap;
    }

    /// 
Solve a QAP to optimality.
    /// 
    public void Solve() {

      _results = new BranchAndBoundResults();
      Stack stack = new Stack();
      stack.Push(new BranchAndBoundNode(_qap));

      while (stack.Count > 0) {
        BranchAndBoundNode node = stack.Pop();

        if (node.LowerBound <= _results.Objective) {
          if (node_is_easy_to_solve) {
            EnumerateSolutions(node, _results);
          }
          else {
            node.LowerBound = GLB.Bound(node.Qap, U);

            if (node.LowerBound <= _results.Objective) {
              foreach (BranchAndBoundNode subNode in Branch(node, U).OrderBy(n => n.LowerBound)) {
                stack.Push(subNode);
              }
            }
          }
        }
      }
    }
  }

  /// 
A subproblem in a branch-and-bound algorithm.
  /// 
  public class BranchAndBoundNode {
    private readonly Qap _qap;
    /// 
The size of the subproblem.
    /// 
    public int Size {
      get { return _qap.Size; }
    }

    /// 
The lower bound of the subproblem.
    /// 
    public double LowerBound { get; set; }

    /// 
The QAP subproblem.
    /// 
    public Qap Qap {
      get { return _qap; }
    }
  }

  /// 
Branch-and-bound solution information.
  /// 
  public class BranchAndBoundResults {
    public int[] p { get; set; }
    public double Objective { get; set; }
  }

Let's work through the pseudocode.  Our stack will contain the current set of subproblems that need to be solved.  We initialize the stack with our QAP instance, wrapping it inside a BranchAndBoundNode data  structure.  We'll develop this data structure as we go, but from the code we can already see two uses:

• It stores the subproblem (a QAP),
• It stores the lower bound computed by the code from my last post.

The other thing we do is allocate a data structure which stores the results: the best assignment found so far, and the best objective value.  In the main loop of Solve() we repeatedly pop nodes from the stack.  We first compare the node's lower bound to the best objective (also called the incumbent).  If the lower bound exceeds the incumbent value, there is no way this subproblem can lead to a better result, so we need not process it.  The incumbent value may have changed since the subproblem was created.  Then we check to see if the subproblem is "easy" - if so, we can simply enumerate over all possible solutions for the subproblem and update the results if needed.  Otherwise, we use our GLB code to compute a bound.  Again, we check the lower bound to see if we can eliminate (or fathom) the subproblem.  If not, we need to subdivide the node into subproblems, and push them onto the stack.

The top-level algorithm itself is not all that complicated!  It's also quite general.  We have already described how to compute bounds, so all we need to do is specify:

• How to enumerate solutions for small problems.
• How to branch - i.e. how to subdivide BranchAndBoundNodes.

The motivation and code for enumerating and branching will be the subject of the next two posts.

Yesterday I introduced a simple OML model for computing project schedules.  The model has the following sections:

• Parameters: input data that comes from Excel.
• Decisions: the values that I want to solve for.
• Constraints: the relationships that decisions and parameters need to satisfy.
• Goals: the objective - what I want to optimize.

The exciting cliffhanger from my previous post was to model precedence links.  Let's think about how I could specify the links in Excel.  One way to do it is to have a table that simply lists the "predecessor" and "successor" tasks.  For example:

Pred	Succ
0	1
1	2

The first row says that task 0 is the predecessor of task 1, i.e. task 1 should start after task 0 finishes.  If we want to bind this data to our model, we should introduce a Set and a Parameter, just like we did for Tasks:

  Parameters[Sets, Links],
  Parameters[Integers, Pred[Links], Succ[Links]],

The domain for Pred and Succ is Integers because the table stores indexes for the tasks.  Now I need to change the model to enforce the links.  That's a constraint, which is expressed straightforwardly using a Foreach:

    Foreach[{l, Links}, Start[Succ[l]] >=Finish[Pred[l]] ]

If I setup the data binding for Links and re-solve, then I will see that the links are honored.  If you compare to Project's schedule, you may spot a problem: some of the tasks (for example t3) are scheduled later than they are in Project.  This is totally valid according to our model - t3 has "slack" and can move around without affecting the project schedule.  Think about how you could fix up the model to match Project, so that all tasks are scheduled "as soon as possible".

On to the "extra credit" problem: how to model "start no earlier than" (SNET) constraints.  This is actually not too hard if we just consider the SNET date to be another parameter just like Duration.  The constraint is straightforward:

  Parameters[Reals, SNET[Tasks]],
  Constraints[ Foreach[{t, Tasks}, Start[t] >= SNET[t] ]]

In my Excel spreadsheet, I added a SNET column with the date for the constraint for each task.  If a task does not have a SNET constraint, I just fill in a dummy date (e.g. 1/1/2009).  But I need a number, not a date, so in the adjacent column I add the formula "=C2-TODAY()" to give me the number of days since the project start date.  That's what I want to bind to.

For reference, here is the data in my spreadsheet:

 

And here are the data bindings I created using the "Binding" button in the Solver Foundation ribbon:

 

If you've been working with the spreadsheet, you may have been annoyed by a couple of things:

• If I add new tasks, I need to change the data binding to include the new cells.
• The conversion between dates and numbers is kind of bothersome.
• This is getting complicated.  For example, Project has eight different constraint types.  If I wanted to model them all the same way I modeled SNET, I would need to add a whole bunch of spreadsheet columns. 

Maybe it's time to turn this OML model into a C# program that directly calls Solver Foundation Services to model and solve project schedules.  Next time we will do just that.

I'd like to work through a few examples of how to model real problems using Solver Foundation.  If you download the Solver Foundation Express Edition, you can follow along: start Excel, click the Solver Foundation tab, then click on the Model button to get a blank modeling pane.

I spent my first five years at Microsoft on the Project team, writing the project scheduling engine for Project Server.  Project is a great team and an awesome product, and I still have a soft spot for project scheduling.  My intent is not to reinvent Project - just to make a connection with a familiar problem so I can show off the Solver Foundation platform while introducing some key concepts. For those who don't know Project, the goal is to come up with a schedule for a list of tasks in a project.  Here are a few of the basics:

• Tasks have a start, finish, and duration.
• No task starts before a given project start date.
• The finish date is computed as start + duration.
• A pair of tasks can be linked by precedence relationships.  When there is a link from a task A to a task B, the B's start must occur on or after A's finish date.
• The project finish date is defined as the latest finish date of any task in the project.

Given these requirements, we would like to compute the start and finish dates for all of the tasks in the project, and we want to find a schedule ends as soon as possible.  Let's use the Solver Foundation to create a model in OML (Solver Foundation's declarative language) and compute some schedules.  I always find it helpful to work out an example before I try and model a problem.  Here's a simple schedule in Project with five tasks, and two links that connect t0, t1, t2.  To compute the schedule by hand I could use a simple critical path scheduling algorithm.

 

Right away you can spot a few issues:

• I haven't taken calendars into account - normally I don't want to schedule work on Saturdays and Sundays.
• I haven't taken resources into account - in real life, people (or machines, or rooms) are assigned to work on tasks, and their availability may affect my schedule.
• There are lots of cool Project features that I haven't taken into account - e.g. task constraints, summary tasks, etc.
• I haven't taken uncertainty into account - in real life, durations aren't known precisely.

I will revisit some of these points in future posts - these are all things that I can address with Solver Foundation.  For now, I am going to stick with a simple version of the problem.  I'd rather get a correct model on a simplified version of the problem than an incorrect model on a more complex one!  In this sense, modeling is no different from writing an essay or designing software.

There are a few questions to ask yourself when building an OML model:

• What are the key objects or entities that I will be working with? (tasks, links)
• What do I want the solver to compute? (start and finish dates for tasks)
• What do I know?  (the durations on the tasks)
• What are my restrictions? (links enforce an ordering on tasks, start + duration = finish)
• What do I want to optimize? (the finish date on the project)

If you have good answers to these questions, and start simple, you can usually come up with an OML model that describes the problem.  Let's start by modeling the start and finish dates of the tasks, since that is ultimately what we want to find out.  Each task has a start and a finish - these are dates.  Solver Foundation does not directly support date values, so let's say that start and finish are real numbers, representing the number of days since the project start date.  To express the idea that each start and finish is associated with a task, we define a "set" Parameter.  The start and finish dates are the values that we want the solver to figure out, so these are Decisions:

Model[
 Parameters[Sets, Tasks],  
 
 Decisions[
     Reals[0, Infinity], Start[Tasks], Finish[Tasks]
  ],

Notice the syntax for Decisions: the first argument indicates the possible values (domain) for the decisions.  We could have used Integers[0, Infinity] for the first argument, but then we would not be able to model fractional days. Loosely speaking, the start and finish dates are the "output".  What's input?  The durations.  Let's suppose that the durations are entered into our spreadsheet like this:

Task	Dur
t0	1
t1	2
t2	3
t3	2
t4	4

We will define Duration as a Parameter so that we can use data binding to determine the values.  Just like Start and Finish, Duration is indexed by Task. 

Parameters[Reals, Duration[Tasks]],

After entering the duration data into Excel I can use the Binding button to add a binding from the range B2:B6 to Duration.  When I solve the model, Solver Foundation will use whatever values are entered in the spreadsheet.  But let's not get ahead of ourselves...next let's express the relation between Start, Finish, and Duration.  Relationships that need to be enforced are described using Constraints.  OML provides a bunch of operators and expressions for describing constraints.  In this case, we have one constraint per task.  This is easy to express using Foreach:

 Constraints[
    Foreach[{t, Tasks}, Start[t] + Duration[t] == Finish[t] ]
 ] 

We almost have a working model - the last thing we need to express is our Goal: that we want the project to finish as soon as possible.  So let's introduce ProjectFinish as another decision:

 Decisions[
     Reals[0, Infinity], Start[Tasks], Finish[Tasks], ProjectFinish
  ]

The ProjectFinish is the maximum finish date of any task in the model.  In other words, it is at least as big as the Finish date of each task, so we add this constraint.

  Foreach[{t, Tasks}, ProjectFinish >= Finish[t] ]

Now we can express our goal: minimize the project finish:

  Goals[ Minimize[finish -> ProjectFinish ] ]

The "->" syntax simply gives a label to our goal, it's optional.  Here's the complete model:

Model[
  Parameters[Sets, Tasks],  
  Parameters[Reals, Duration[Tasks]],

  Decisions[
     Reals[0, Infinity], Start[Tasks], Finish[Tasks], ProjectFinish
  ],

  Constraints[
    Foreach[{t, Tasks}, Start[t] + Duration[t] == Finish[t] ],
    Foreach[{t, Tasks}, ProjectFinish >= Finish[t] ],
  ],
 
  Goals[ Minimize[ finish -> ProjectFinish ]  ]
]

If you have entered this into the modeling pane, and used the Binding button to define the data binding between the spreadsheet data and Duration, you can hit Solve and observe the Start and Finish dates.  You will see something like this:

 

If you want to turn the starts and finishes into dates, you can enter a formula like "=now() + C5" in the adjacent column.  Actually I put this formula in a column next to my duration values, so I can see the durations, starts, and finishes all in one place.

You should be thoroughly unimpressed by this model (as I am ;)), because it doesn't do all that much.  All I'm doing is letting Solver Foundation compute finish = start + duration, which I could have done very easily in Excel.  I need to model "links": a link causes the start date of the successor task to be later than the finish of the predecessor task.  If you've made it this far, you have all the building blocks you need to model this in OML.  Feel free to post your solution here, or just wait until tomorrow for my solution.

Extra credit: if you use Project, you know that you can define "constraints" on tasks.  For example, I can double click on a task and say that a task must "Start No Earlier Than March 30, 2009".  These constraints need to be honored, in addition to links.  How can I model SNET constraints using OML?

 

Whew - an exciting and tiring week is over.  We had hundreds of employees stop by our TechFest booth to ask questions, talk about their team's optimization or modeling problems, look at our demos, or simply hear more about Solver Foundation.   A big thanks to my teammates Lengning, Lucas, and Lin for doing such a great job this last week. 

OML - our declarative language for specifying optimization problems - was received with great interest.  Microsoft France's experience using Solver Foundation to schedule TechDays really "worked" as an example. Even our simple quadratic programming example in Excel opened a lot of eyes.  There several questions about data binding (and there have been some recent threads in our forums), so I think I'd like to devote a couple of posts to that.  Lastly, I was surprised to get several inquiries into our nonlinear unconstrained optimizer - that's an area I will devote some time to that as well.

I will be attending a couple of conferences over the next several months, but for now it's back to work on the next version of Solver Foundation.

Another busy day - we spoke with literally hundreds of Microsoft employees about Solver Foundation at TechFest '09.  Microsoft is a huge company, so we got a chance to talk to people from many different parts of the company with varying backgrounds.  Some are intimately familiar with optimization, others are not.  So part of the challenge is to try to "meet people where they are at" and talk about Solver Foundation in a way that makes sense for them, without trivializing it or making it sound like magic beans.  If you had the opportunity to come by our booth tomorrow, you'd get a few cool demos, and a variation on the following "pitch".  Depending on what you wanted to talk about we could talk more about how to model real problems, or maybe the API, or maybe some of the solvers that we're developing.  Or maybe the stock market, or how it sucked to lose the Sonics.  If people leave with a good idea of what we're about and can think of some situations where Solver Foundation could help, I'm happy.  Anyway, here's goes:

People and businesses need to be able to juggle different priorities and constraints to make good decisions.  Examples include production planning, scheduling projects, configuring IT systems, advertising online.  Microsoft Solver Foundation is a managed code platform for planning, scheduling, configuration, and optimization.   Solver Foundation lets you easily describe your problem, get it solved, and connect it to your application. 

It has three main layers:
 • Modeling and Programming: the modeler lets you describe optimization problems declaratively: "what", not "how".  You can do it in code, or in our Excel add-in.
 • Solver Foundation Services: a full-featured .Net library that can be accessed from Visual Studio, C#, VB, ASP.Net, Silverlight.  It transparently handles parallelism & multiple cores.  It provides a rich object model, events, and data binding.
 • Solvers: we provide a bunch of solvers that cover a wide range of optimization problems including: LP, QP, constraint, MIP, nonlinear unconstrained.  We feature an open, extensible architecture that lets plug-in third party solvers if you so choose.

Solver Foundation is developed in conjunction with researchers in Redmond and Cambridge, and is a great example of Microsoft Research innovation.  We recently released our 1.1 version and you can check out solverfoundation.com to download our free Express version!