// (c) 2023-2024 Fair Isaac Corporation
using System;
using System.Linq;
using Optimizer;
using Optimizer.Objects;
using static Optimizer.Objects.Utils;
using static Optimizer.Objects.ConstantExpression;
namespace XpressExamples
{
/// <summary>Cutting stock problem</summary>
/// <remarks>
/// Solved by column (= cuttingpattern) generation heuristic looping
/// over the root node.
/// </remarks>
internal class CuttingStock
{
const int NWIDTHS = 5;
const int MAXWIDTH = 94;
const double EPS = 1e-6;
const int MAXCOL = 10;
/****Data****/
static readonly double[] width = { 17, 21, 22.5, 24, 29.5 }; /* Possible widths */
static readonly int[] demand = { 150, 96, 48, 108, 227 }; /* Demand per width */
/****Shared Variables****/
static Variable[] vPattern; /* Rolls per pattern */
static Inequality[] cDemand; /* Demand constraints */
static LinTermMap eObj; /* Objective function */
private static void PrintProblemSolution(XpressProblem prob)
{
double[] sol = prob.GetSolution();
Console.WriteLine("Objective: " + prob.ObjVal);
Console.WriteLine("Variables:");
Console.Write("\t");
foreach (Variable v in prob.GetVariables())
{
Console.Write("[{0}: {1:f2}] ", v.GetName(), v.GetValue(sol));
}
Console.WriteLine();
}
/**************************************************************************/
/* Integer Knapsack Algorithm for solving the integer knapsack problem */
/* z = max{cx : ax <= R, x <= d, x in Z^N} */
/* */
/* Input data: */
/* N: Number of item types */
/* c[i]: Unit profit of item type i, i=1..n */
/* a[i]: Unit resource use of item type i, i=1..n */
/* R: Total resource available */
/* d[i]: Demand for item type i, i=1..n */
/* Return values: */
/* xbest[i]: Number of items of type i used in optimal solution */
/* zbest: Value of optimal solution */
/**************************************************************************/
private static double SolveKnapsack(int N, double[] c, double[] a, double R, int[] d, int[] xbest)
{
using (XpressProblem prob = new XpressProblem())
{
Variable[] x = prob.AddVariables(N)
.WithType(ColumnType.Integer)
.WithName(i => $"x_{i}")
.WithLB(i => d[i])
.ToArray();
prob.AddConstraint(
ScalarProduct(x, a) <= R
);
prob.SetObjective(
ScalarProduct(x, c),
Optimizer.ObjSense.Maximize
);
prob.MipOptimize();
double zbest = prob.ObjVal;
double[] sol = prob.GetSolution();
xbest = sol.Select(value => (int)Math.Round(value)).ToArray();
return zbest;
}
}
private static void ModCutStock(XpressProblem p)
{
// (Basic) cutting patterns
int[,] patterns = new int[NWIDTHS, NWIDTHS];
foreach (int j in Enumerable.Range(0, NWIDTHS))
{
patterns[j, j] = (int)Math.Floor((double)MAXWIDTH / width[j]);
}
vPattern = p.AddVariables(NWIDTHS)
.WithUB(i => Math.Ceiling((double)demand[i] / patterns[i, i]))
.WithName(i => $"pat_{i + 1}")
.ToArray();
// Minimize total number of rolls
p.SetObjective(Sum(vPattern));
// Satisfy the demand per width;
// Sum of patterns must exceed demand
cDemand = p.AddConstraints(NWIDTHS,
i => (Sum(NWIDTHS, j => vPattern[j] * patterns[i, j]) >= demand[i]).SetName($"Demand_{i}")).ToArray();
}
/**************************************************************************/
/* Column generation loop at the top node: */
/* solve the LP and save the basis */
/* get the solution values */
/* generate new column(s) (=cutting pattern) */
/* load the modified problem and load the saved basis */
/**************************************************************************/
private static void SolveCutStock(XpressProblem p)
{
// Basis information
int[] rowstat = new int[p.Rows];
int[] colstat = new int[p.Cols];
int[] x = new int[NWIDTHS];
long starttime = Environment.TickCount64;
foreach (int npass in Enumerable.Range(0, MAXCOL))
{
p.LpOptimize(); // Solve the LP
if (p.SolStatus != Optimizer.SolStatus.Optimal)
throw new Exception("failed to optimize with status " + p.SolStatus);
p.GetBasis(rowstat, colstat); // Save the current basis
// solve integer knpsack problem
double z = SolveKnapsack(NWIDTHS, p.GetDuals(cDemand), width, MAXWIDTH, demand, x);
Console.WriteLine("Iteration {0:d}, {1:f2} sec",
npass + 1, (Environment.TickCount64 - starttime) / 1000.0
);
if (z < 1 + EPS)
{
Console.WriteLine("No profitable column found.");
Console.WriteLine();
}
else
{
Console.WriteLine("New pattern found with marginal cost {0}", z - 1);
// Create a new variable for this pattern:
Variable v = p.AddVariable(
0.0,
Double.PositiveInfinity,
ColumnType.Integer,
String.Format("pat_{0:d}", p.Cols + 1)
);
// Add new variable to the objective
v.ChgObj(1.0);
// Add new variable to demand constraints
double ub = 0.0;
for (int i = 0; i < NWIDTHS; ++i)
{
if (x[i] > 0)
{
cDemand[i].ChgCoef(v, x[i]);
ub = Math.Max(ub, Math.Ceiling(demand[i] / (double)x[i]));
}
}
// Change the upper bound on the new variable
v.SetUB(ub);
p.LoadBasis(rowstat, colstat);
}
}
p.MipOptimize();
if (p.SolStatus != Optimizer.SolStatus.Optimal)
throw new Exception("failed to optimize with status " + p.SolStatus);
PrintProblemSolution(p);
}
public static void Main(string[] args)
{
using (XpressProblem prob = new XpressProblem())
{
// Output all messages.
prob.callbacks.AddMessageCallback(DefaultMessageListener.Console);
prob.OutputLog = 0;
prob.MIPLog = 0;
ModCutStock(prob);
SolveCutStock(prob);
}
}
}
}
|
