Cutting stock: column generation implemented by C++

这篇文章是讲述如何用C++来实施列生成算法并求解钢卷切割问题. 这篇文章来自Cplex的用户手册中Discrete Optimization的Cutting Stock: Column Generation. 为保持原意，我先贴原文，然后再翻译解读. 由于我也是初学，因此难免会意识不到自己的翻译错误。如果读者发现了，please don’t hesitate to tell me. Thank you! 废话少说开始正题
This chapter uses an example of cutting stock to demonstrate the technique of column generation in Concert Technology. In it, you will learn:
◆ how to use classes of ILOG CPLEX for column generation in column-wise modeling;
◆ how to modify a model and re-optimize;
◆ how to change the type of a variable with IloConversion;
◆ how to use more than one model;
◆ how to use more than one algorithm (instances of IloCplex, for example).
在本文将会学到以下五种技术
（1）再列生成模型中，如何使用cplex 的类生成列；
（2）如何改变一个模型并重新优化该模型；
（3）如何使用 IloConversion 改变一个变量的类型；
（4）在一个算法中如何使用多个模型；
（5）如何使用调用多个算法。

What Is Column Generation?

In colloquial terms, column generation is a way of beginning with a small, manageable part of a problem (specifically, a few of the variables), solving that part, analyzing that partial solution to determine the next part of the problem (specifically, one or more variables) to add to the model, and then resolving the enlarged model. Column-wise modeling repeats that process until it achieves a satisfactory solution to the whole of the problem.
总的来说，列生成算法是从主问题的一部分（子问题）出发，子问题是一个小的，可调整的问题。然后，求解子问题，分析该问题的解，然后找到下一个子问题的解，并添加到主问题的模型中，不断求解，直到产生一个主问题的满意解。
解读：列生成算法的思想是分治算法，把一个运算规模很大的问题，分解成很多个可解的子问题。并且这些的子问题无需全部列出来。先从一个可行的子问题入手，求得一个解。然后把解返回到主问题中，主问题根据反馈的解再产生子问题。这个过程不断重复，直到主问题的最优解（可接受的解）能够求得。

In formal terms, column generation is a way of solving a linear programming problem that adds columns (corresponding to constrained variables) during the pricing phase of the simplex method of solving the problem. In gross terms, generating a column in the primal simplex formulation of a linear programming problem corresponds to adding a constraint in its dual formulation. In the dual formulation of a given linear programming problem, you might think of column generation as a cutting plane method.

In that context, many researchers have observed that column generation is a very powerful technique for solving a wide range of industrial problems to optimality or to near optimality. Ford and Fulkerson, for example, suggested column generation in the context of a multi-commodity network flow problem as early as 1958 in the journal of Management Science. By 1960, Dantzig and Wolfe had adapted it to linear programming problems with a decomposable structure. Gilmore and Gomory then demonstrated its effectiveness in a cutting stock problem. More recently, vehicle routing, crew scheduling, and other
integer-constrained problems have motivated further research into column generation.

Column generation rests on the fact that in the simplex method, the solver does not need access to all the variables of the problem simultaneously. In fact, a solver can begin work with only the basis (a particular subset of the constrained variables) and then use reduced
cost to determine which other variables to access as needed.

Column-Wise Models in Concert Technology

Concert Technology offers facilities for exploiting column generation. In particular, you can design the model of your problem (one or more instances of the class IloModel) in terms of a basis and added columns (instances of IloNumVar, IloNumVarArray, IloNumColumn, or IloNumColumnArray). For example, instances of IloNumColumn represent columns,and you can use operator() in the classes IloObjective and IloRange to create terms in column expressions. In practice, the column serves as a kind of place holder for a variable in other extractable objects (such as a range constraint or an objective) when your application needs to declare or use those other extractable objects before it can actually know the value of a variable appearing in them.

Furthermore, an instance of IloCplex provides a way to solve the master linear problem, while other Concert Technology algorithms (that is, instances of IloSolver, of IloCplex itself, or of subclasses of IloAlgorithm,for example) lend themselves to other parts of the problem by determining which variables to consider next (and thus which columns to generate).
In the Concert Technology Reference Manual, the concept Column-Wise Modeling provides more detail about this topic and offers simple examples of its use.

Describing the Problem

The cutting stock problem in this chapter is sometimes known in math programming terms as a knapsack problem with reduced cost in the objective function.

Generally, a cutting stock problem begins with a supply of rolls of material (the stock). Strips are cut from a continuous roll. All the strips cut from one roll are known together as a pattern. The point of this example is to use as few rolls of stock as possible to satisfy some specified demand. By convention, it is assumed that only one pattern is laid out across the stock; consequently, only one dimension—the width—of each roll of stock is important.
Even with that simplifying assumption, the fact that there can be so many different patterns makes a naive model of this problem (where a user declares one variable for every possible pattern) impractical. Such a model introduces too many symmetries. Fortunately, for any given customer order, a limited number of patterns will suffice, so many of the possible
patterns can be disregarded, and the application can focus on finding the relevant ones. Here is a conventional statement of a cutting stock problem in terms of the unknown Xj, the number of times that pattern j will be used, and Aij, the number of items i of each pattern j needed to satisfy demand di:

Master Problem:

Minimize Sum(j) : X[j];
Subject to: Sum(i, j): A[i][j] * X[j]>= d[i], for any iterm i;

Solving this model with all columns present from the beginning is practically impossible. In fact, even with only 10 types of items with a size roughly 1/10 of the width of the roll, there would exist roughly 10^10 kinds of patterns, and hence that many decision variables. Such a formulation might not even fit in memory on a reasonably large computer. Moreover, most of those patterns would obviously not be interesting in a solution. These considerations make column generation an interesting pproach for this problem.
To solve a cutting stock problem by column generation, start with a subproblem. Choose one pattern, lay it out on the stock, and cut as many items as possible, subject to the constraints of demand for that item and the width of the stock. This procedure will surely work in that it produces some answer (a feasible solution) to the problem, but it will not necessarily
produce a satisfactory answer in this way.
To move closer to a satisfactory solution, the application can then generate other columns. That is, other decision variables (other Xj) will be chosen to add to the model. Those decision variables are chosen on the basis of their favorable reduced cost. In linear programming terms, choose variables for which the negative reduced cost improves the minimization. In this formulation of the rest of the problem, πi represents the dual value of constraint i, like this:
Sub Problem:

Minimize 1 - sum(i):  π[i]* A[i]
Subject to： sum(i): W[i] * A[i] <= W;

where W is the width of a roll of stock and Ai must be a nonnegative integer. A is a pattern which consists of A[i]. So A[i] is a decision varible.
(The author’s words:The economic meaning of πi is that how many pattern j are needed if the demand of iterm i increases 1; obviously, πi is less than 1/A(i))

Representing the Data

As usual in a Concert Technology application, an environment, an instance of IloEnv, is created first to organize the data and build the model of the problem. The data defining this problem includes the width of a roll of stock.This value is read from a file and represented by a numeric value, rollWidth （scale value W）. The widths of the ordered rolls are also read from a file and put into an array of numeric values, size (array A). Finally, the number of rolls ordered of each width is read from a file and put into an array of numeric values, amount (array d).

Developing the Model: Building and Modifying

In this problem, an initial model cutOpt（the master problem） is built first. Later, through its modifications, another model patGen (the sub problem) is built to generate the new columns.
The first model cutOpt, an instance of IloModel, is declared like this:

IloModel cutOpt (env);

An array of numeric variables (an instance of IloNumVarArray) Cut[j] (X(j)) is declared to represent the decision variables in the model. The index j indicates a pattern using the widths found in the array size (A(i,j)). Consequently, a given variable Cut[j] in the solution will be the number of rolls to cut according to pattern j.

As a model for this problem is built, there will be opportunities to demonstrate to you how to modify a model by adding extractable objects, adding columns, changing coefficients in an objective function, and changing the type of a variable. When you modify a model by means of (by means of = using) the methods of extractable objects, Concert Technology notifies the algorithms (instances of subclasses of IloAlgorithm, such as IloCplex or IloSolver) about the modification. (For more about that idea, see the concept of Notification in the Concert Technology Reference Manual.)

When IloCplex, for example, is notified about a change in an extractable object that it has extracted, it maintains as much of the current solution information as it can accurately and reasonably. Other parts of the ILOG CPLEX User’s Manual offer more detail about how the algorithm responds to modifications in the model.

Adding Extractable Objects: Both Ways

In a Concert Technology application, there are two ways of adding extractable objects to a model: by means of a template function (IloAdd) or by means of a method of the model (IloModel::add or IloArray::add). In this example, you see both ways.

Using a Template to Add Objects

When an objective is added to the model, the application needs to keep a handle to the objective RollsUsed because it is needed when the application generates columns. For that purpose, the application relies on the template function IloAdd, like this:

IloObjective RollsUsed = IloAdd(cutOpt, IloMinimize(env));

Apart from the fact that it preserves type information, that single line is equivalent to these lines:

IloObjective RollsUsed = IloMinimize(env);
cutOpt.add(RollsUsed);

Likewise, the application adds an array of constraints to the model. These constraints are needed later in column generation as well, so the application again uses IloAdd again to add the array Fill to the model.

IloRangeArray Fill = IloAdd(cutOpt, IloRangeArray(env, amount, IloInfinity));

That statement creates amount.getSize range constraints. Constraint Fill[i] has a lower bound of amount[i] and an upper bound of IloInfinity
We can see that Fill is the constraint array.
Fill[i] = Sum(i, j): A[i][j] * X[j]>= d(i);

Using a Method to Add Objects

It is also possible to add objects to your model by means of the method IloModel::add. This example uses that approach in these lines:

patGen.add(IloScalProd(size, Use) <= rollWidth);

The author: This line is the constraint set. Where is the Use?
In contrast, to create Cut, an array of decision variables, the application uses the add method of IloArray within a loop, like this:

IloInt nWdth = size.getSize();
for (j = 0; j < nWdth; j++)
     Cut.add( IloNumVar(  RollsUsed(1) + Fill[ j ]( int(rollWidth / size[ j ]) )  )  );

In that loop, the variable Cut[j] (X(j)) is created with an objective coefficient of 1 (one), and all its pattern coefficients are initially 0 (zero) except for the one in Fill[j]. That one is set to (int(rollWidth / size[j])). By default, the type of each variable is IloFloat, so it is not explicitly mentioned in the loop.

The author’s words:
RollsUsed is the objective function, Fill[ j ] is the constraint [ j ], rollWidth is W, size is A[i, j], int (); Each Cut only generate a type of size of parttern. It is a initial and feasible solution. For example, for a roll with width 10, an three type items are needed, 3, 5, 7. size. getSize() calculates the types of items.Therefore, in the initial solution. a roll is cut into three items with width 3 , or two items with width 5, or a item with width 7.

Adding Columns to a Model

Creating a new column to add to a model in Concert Technology is a two-step process:

1. Create a column expression defining the new column.
2. Create a variable using that column expression and add the variable to the model.
   For example, in this problem, RollsUsed is an instance of IloObjective. The statement RollsUsed(1) creates a term in a column expression defining how to add a new variable as a linear term with a coefficient of 1 (one) to the expression RollsUsed.
   The terms of a column expression are connected to one another by the overloaded operator +.The constraints of this problem are represented in the array Fill. A new column to be added to the model has a coefficient for each of the constraints in Fill. Those coefficients are represented in the array newPatt.

Changing Coefficients of an Objective Function

According to that two-step procedure for adding a column to a model, the following lines create the column with coefficient 1 (one) for the objective RollsUsed and create the new variable with bounds at 0 (zero) and at MAXCUT.

IloNumColumn col = RollsUsed(1);
for (IloInt i = 0; i < Fill.getSize(); ++i)
col += Fill[i](newPatt[i]);
IloNumVar var(col, 0, MAXCUT);

(However, those lines do not appear in the example at hand.) Concert Technology offers a shortcut in the operator() for an array of range constraints. Those lines of code can be condensed into the following line:

IloNumVar var(RollsUsed(1) + Fill(newPatt), 0, MAXCUT);

In other words, Fill(newPatt) returns the column expression that the loop would create. You will see a similar shortcut in the example.
There are also other ways of modifying an objective in a Concert Technology model. For example, you can call the method IloObjective::setCoef to change a coefficient in an objective. In this example, you see this way of changing a coefficient when the application wants to find a new pattern.

ReducedCost.setCoef(Use, price);

You can also change the sense of an objective in Concert Technology. There are at least two ways to do so. One approach uses the method setSense of the class IloObjective. The other approach first removes the current objective function from the model, then uses either the function IloMinimize or IloMaximize to create a new objective function, and finally adds that new objective to the model.

Changing the Type of a Variable

With Concert Technology, in order to change the type of a variable in a model, you actually create an extractable object (an instance of IloConversion, documented in the Concert Technology Reference Manual) and add that object to the model.
In the example, when the application needs to change the elements of Cut (an array of numeric variables) from their default type of ILOFLOAT to integer (type ILOINT), it creates an instance of IloConversion, applies it the array Cut, and adds the conversion to the model, cutOpt, like this:

cutOpt.add(IloConversion(env, Cut, ILOINT));

Cut Optimization Model

Here is a summary of the initial model cutOpt:

IloModel cutOpt (env);
IloObjective RollsUsed = IloAdd(cutOpt, IloMinimize(env));
IloRangeArray Fill = IloAdd(cutOpt, IloRangeArray(env, amount, IloInfinity));
IloNumVarArray Cut(env);
IloInt nWdth = size.getSize();
for (j = 0; j < nWdth; j++)
   Cut.add(IloNumVar(RollsUsed(1) + Fill[j](int(rollWidth / size[j]))));

Pattern Generator Model

The application in this example starts with an initial model cutOpt and through modifications eventually builds a second model patGen. This pattern generator patGen (in contrast to cutOpt) includes integer variables in the array Use. That array appears in the only constraint added to patGen: a scalar product making sure that the patterns used do not exceed the width of rolls. The application also adds a rudimentary objective function to patGen. This objective initially consists of only the constant 1 (one). The rest of the objective function depends on the solution found with the initial model cutOpt. The application will build that objective function as that information is computed. Here, in short, is patGen:

IloModel patGen (env);
IloObjective ReducedCost = IloAdd(patGen, IloMinimize(env, 1));
IloNumVarArray Use(env, nWdth, 0, IloInfinity, ILOINT);
patGen.add(IloScalProd(size, Use) <= rollWidth);

Solving the Problem: Using More than One Algorithm

This example does not solve the problem to optimality. It only generates a good feasible solution. It does so by first solving a linear relaxation of the problem. In other words, the decision variables Cut[j] are relaxed (not required to be integers). This linear relaxation of the problem then generates columns. While generating columns, the application keeps the variables generated so far, changes their type to integer, and resolves the problem. As you’ve seen, this example effectively builds two models of the problem, cutOpt and patGen. Likewise, it uses two algorithms (that is, two instances of IloCplex) to arrive at a
good feasible solution. Here’s how to create the first algorithm cutSolver and extract the initial model cutOpt:

IloCplex cutSolver(cutOpt);

And here is how to create the second algorithm and extract patGen:

IloCplex patSolver(patGen);

The heart of the example is here, in the column generation and optimization over current patterns:

IloNumArray price(env, nWdth);
IloNumArray newPatt(env, nWdth);
for (; ; ) {
/// OPTIMIZE OVER CURRENT PATTERNS ///
	cutSolver.solve();
	report1 (cutSolver, Cut, Fill);
/// FIND AND ADD A NEW PATTERN ///
	for (i = 0; i < nWdth; i++)
		price[i] = -cutSolver.getDual(Fill[i]);
	ReducedCost.setCoef(Use, price);
	patSolver.solve();
	report2 (patSolver, Use, ReducedCost);
	if (patSolver.getValue(ReducedCost) > -RC_EPS) break;
	patSolver.getValues(newPatt, Use);
	Cut.add( IloNumVar(RollsUsed(1) + Fill(newPatt)) );
}
cutOpt.add(IloConversion(env, Cut, ILOINT));
cutSolver.solve();

Those lines solve the current subproblem cutOpt by calling cutSolver.solve. Then they copy the values of the negative dual solution into the array price. They use that array to set objective coefficients in the model patGen. Then they solve the right pattern generation problem.
If the objective value exceeds the negation of the optimality tolerance (-RC_EPS), then the application has proved that the current solution of the model cutOpt is optimal within the given optimality tolerance (RC_EPS). Otherwise, the application copies the solution of the current pattern generation problem into the array newPatt and uses that new pattern to build
the next column to add to the model cutOpt. Then it repeats the procedure.

Ending the Program

As in other C++ Concert Technology applications, this program ends with a call to IloEnv::end to de-allocate the models and algorithms once they are no longer in use.

env.end();

Complete Program

You can see the entire program here or online in the standard distribution of ILOG CPLEX in the file /examples/src/cutstock.cpp. To run the example, you need only a license for ILOG CPLEX. It is also possible to write your own pattern generation models to use with algorithms of ILOG Solver and ILOG Hybrid Optimizers. In that case, you need a
license for those products.

// -------------------------------------------------------------- -*- C++ -*-
// File: examples/src/cutstock.cpp
// Version 9.0
// --------------------------------------------------------------------------
// Copyright (C) 1999-2003 by ILOG.
// All Rights Reserved.
// Permission is expressly granted to use this example in the
// course of developing applications that use ILOG products.
// --------------------------------------------------------------------------
//
#include <ilcplex/ilocplex.h>
ILOSTLBEGIN
#define RC_EPS 1.0e-6
typedef IloArray<IloNumArray> IloNumArray2;
static void readData (const char* filename, IloNum& rollWidth,
IloNumArray& size, IloNumArray& amount);
static void report1 (IloCplex& cutSolver, IloNumVarArray Cut, IloRangeArray Fill);
static void report2 (IloAlgorithm& patSolver,
IloNumVarArray Use,
IloObjective obj);
static void report3 (IloCplex& cutSolver, IloNumVarArray Cut);
/// MAIN PROGRAM ///
int
main(int argc, char * *argv)
{
	IloEnv env;
	try {
		IloInt i, j;
		IloNum rollWidth;
		IloNumArray amount(env);
		IloNumArray size(env);
		if ( argc > 1 )
			readData(argv[1], rollWidth, size, amount);
		else
			readData("../../../examples/data/cutstock.dat", rollWidth, size, amount);
/// CUTTING-OPTIMIZATION PROBLEM ///
		IloModel cutOpt (env);
		IloObjective RollsUsed = IloAdd(cutOpt, IloMinimize(env));
		IloRangeArray Fill = IloAdd(cutOpt, IloRangeArray(env, amount, IloInfinity));
		IloNumVarArray Cut(env);
		IloInt nWdth = size.getSize();
		for (j = 0; j < nWdth; j++) {
			Cut.add(IloNumVar(RollsUsed(1) + Fill[j](int(rollWidth / size[j]))));
		}
		IloCplex cutSolver(cutOpt);
		/// PATTERN-GENERATION PROBLEM ///
		IloModel patGen (env);
		IloObjective ReducedCost = IloAdd(patGen, IloMinimize(env, 1));
		IloNumVarArray Use(env, nWdth, 0, IloInfinity, ILOINT);
		patGen.add(IloScalProd(size, Use) <= rollWidth);
		IloCplex patSolver(patGen);
/// COLUMN-GENERATION PROCEDURE ///
		IloNumArray price(env, nWdth);
		IloNumArray newPatt(env, nWdth);
/// COLUMN-GENERATION PROCEDURE ///
		for (;;) {
/// OPTIMIZE OVER CURRENT PATTERNS ///
			cutSolver.solve();
			report1 (cutSolver, Cut, Fill);
/// FIND AND ADD A NEW PATTERN ///
			for (i = 0; i < nWdth; i++) {
				price[i] = -cutSolver.getDual(Fill[i]);
			}
			ReducedCost.setCoef(Use, price);
			patSolver.solve();
			report2 (patSolver, Use, ReducedCost);
			if (patSolver.getValue(ReducedCost) > -RC_EPS) break;
			patSolver.getValues(newPatt, Use);
			Cut.add( IloNumVar(RollsUsed(1) + Fill(newPatt)) );
		}
		cutOpt.add(IloConversion(env, Cut, ILOINT));
		cutSolver.solve();
		report3 (cutSolver, Cut);
	}
	catch (IloException& ex) {
		cerr << "Error: " << ex << endl;
	}
	catch (...) {
		cerr << "Error" << endl;
	}
	env.end();
	return 0;
}
static void readData (const char* filename, IloNum& rollWidth,
IloNumArray& size, IloNumArray& amount)
{
	ifstream in(filename);
	if (in) {
		in >> rollWidth;
		in >> size;
		in >> amount;
	}
	else {
		cerr << "No such file: " << filename << endl;
		throw(1);
	}
}
static void report1 (IloCplex& cutSolver, IloNumVarArray Cut, IloRangeArray Fill)
{
	cout << endl;
	cout << "Using " << cutSolver.getObjValue() << " rolls" << endl;
	cout << endl;
	for (IloInt j = 0; j < Cut.getSize(); j++) {
	cout << " Cut" << j << " = " << cutSolver.getValue(Cut[j]) << endl;
	}
	cout << endl;
	for (IloInt i = 0; i < Fill.getSize(); i++) {
		cout << " Fill" << i << " = " << cutSolver.getDual(Fill[i]) << endl;
	}
	cout << endl;
}
static void report2 (IloAlgorithm& patSolver, IloNumVarArray Use, IloObjective obj)
{
	cout << endl;
	cout << "Reduced cost is " << patSolver.getValue(obj) << endl;
	cout << endl;
	if (patSolver.getValue(obj) <= -RC_EPS) {
		for (IloInt i = 0; i < Use.getSize(); i++) {
			cout << " Use" << i << " = " << patSolver.getValue(Use[i]) << endl;
		}
		cout << endl;
	}
}
static void report3 (IloCplex& cutSolver, IloNumVarArray Cut)
{
	cout << endl;
	cout << "Best integer solution uses "<< cutSolver.getObjValue() << " rolls" << endl;
	cout << endl;
	for (IloInt j = 0; j < Cut.getSize(); j++) {
		cout << " Cut" << j << " = " << cutSolver.getValue(Cut[j]) << endl;
	}
}