F# Linear Programming Example

#light

namespace CenterSpace.NMath.Examples.FSharp

open System

open CenterSpace.NMath.Core

  /// <summary>
  /// A .NET example in F# showing how to solve a linear system using linear programming and
  /// the simplex method.
  /// </summary>
module LinearProgramming =
      // A farmer has 640 acres of farmland. It can be planted with wheat, barley, corn or a
      // combination of the three. The farmer wishes to maximize his profit subject to the
      // limits on land, fertilizer, and water.

      // Currently, wheat is $3.38/bushel. The farmer can expect a yield of 55 bushels/acre.
      let wheatPrice = 3.38
      let wheatYield = 55.0
      let wheatRevenuePerAcre = wheatPrice * wheatYield

      // Currently, barley is $1.98/bushel. The farmer can expect a yield of 75 bushels/acre.
      let barleyPrice = 1.98
      let barleyYield = 75.0
      let barleyRevenuePerAcre = barleyPrice * barleyYield

      // Currently, corn is $1.70/bushel. The farmer can expect a yield of 110 bushels/acre.
      let cornPrice = 1.70
      let cornYield = 110.0
      let cornRevenuePerAcre = cornPrice * cornYield

      // Therefore, the objective function is:
      printfn "Maximize"
      printfn "%Aw + %Ab + %Ac" wheatRevenuePerAcre barleyRevenuePerAcre cornRevenuePerAcre
      printfn "where"

      let revenue = new DoubleVector(wheatRevenuePerAcre, barleyRevenuePerAcre, cornRevenuePerAcre)

      // Make a matrix big enough for 5 constraints and 3 variables.
      let constraints = new DoubleMatrix(5, 3)

      // Make a vector of right-hand sides.
      let rightHandSides = new DoubleVector(constraints.Rows)

      // The farmer has 8,000 lbs of nitrogen fertilizer. Its known that wheat requires
      // 12 lb/acre, barley 5 lb/acre and corn 22 lb/acre.
      printfn "12w + 5b + 22c <= 8000"
      constraints.[0, Slice.All] <- new DoubleVector(12.0, 5.0, 22.0)
      rightHandSides.[0] <- 8000.0

      // The farmer has 22,000 lbs of phosphate fertilizer. Its known that wheat requires
      // 30 lb/acre, barley 8 lb/acre and corn 50 lb/acre.
      printfn "30w + 8b + 50c <= 22000"
      constraints.[1, Slice.All] <- new DoubleVector(30.0, 8.0, 50.0)
      rightHandSides.[1] <- 22000.0

      // The farmer has a permit for 1,000 acre-feet of water. Wheat requires 1.5 ft of water, 
      // barley requires 0.7 and corn 2.2.
      printfn "1.5w + 0.7b + 2.2c <= 1200"
      constraints.[2, Slice.All] <- new DoubleVector(1.5, 0.7, 2.2)
      rightHandSides.[2] <- 1200.0

      // The farmer has a maximum of 640 acres for planting.
      printfn "w + b + c <= 640"
      constraints.[3, Slice.All] <- new DoubleVector(1.0, 1.0, 1.0)
      rightHandSides.[3] <- 640.0

      // Create an LP solver with an error tolerance of 0.001.
      let solver = new SimplexLPSolver(0.001)

      // Solve
      solver.Solve(revenue, constraints, rightHandSides, 5, 0, 0) |> ignore

      // Was a finite solution found?
      printfn ""
      if solver.IsGood = true then
        printfn "solution: %s" (solver.Solution.ToString("f0"))
        printfn "optimal value: %s" (solver.OptimalValue.ToString("f0"))
        printfn ""  
      
      // Lets say the farmer is also contractually obligated to farm at least 50 acres of barley.
      printfn "Add variable bound: b >= 10"
      printfn ""  
      let lowerBounds = new DoubleVector(0.0, 10.0, 0.0)
      let upperBounds = new DoubleVector(640.0, 640.0, 640.0)

      // Solve again
      solver.Solve(revenue, constraints, rightHandSides, 5, 0, 0, lowerBounds, upperBounds) |> ignore

      // Good?
      if solver.IsGood = true then
        printfn "solution: %s" (solver.Solution.ToString("f0"))
        printfn "optimal value: %s" (solver.OptimalValue.ToString("f0"))
        printfn ""  

      printfn "Press Enter Key"
      Console.Read() |> ignore