Prev Next Index-> contents reference index search external Up-> CppAD ADFun abs_normal simplex_method simplex_method.hpp ADFun-> record_adfun drivers Forward Reverse sparsity_pattern sparse_derivative optimize abs_normal FunCheck check_for_nan abs_normal-> abs_normal_fun abs_print_mat abs_eval simplex_method lp_box abs_min_linear min_nso_linear qp_interior qp_box abs_min_quad min_nso_quad simplex_method-> simplex_method.cpp simplex_method.hpp simplex_method.hpp Headings

$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$
simplex_method Source Code
namespace CppAD { // BEGIN_CPPAD_NAMESPACE

// BEGIN PROTOTYPE
template <class Vector>
bool simplex_method(
size_t        level   ,
const Vector& A       ,
const Vector& b       ,
const Vector& c       ,
size_t        maxitr  ,
Vector&       xout    )
// END PROTOTYPE
{     // number of equations
size_t ne  = b.size();
// number of x variables
size_t nx = c.size();
CPPAD_ASSERT_UNKNOWN( size_t(A.size()) == ne * nx );
//
if( level > 0 )
{     std::cout << "start simplex_method\n";
}
//
// variables (columns) in the Tableau:
// x: the original primary variables with size n
// s: slack variables, one for each equation
// a: auxillary variables, one for each negative right hand size
// r: right hand size for equations
//
// Determine number of auxillary variables
size_t na = 0;
for(size_t i = 0; i < ne; i++)
{     if( b[i] > 0.0 )
++na;
}
// size of columns in the Tableau
size_t nc = nx + ne + na + 1;

// number of rows in Tableau, the equations plust two objectives
size_t nr = ne + 2;

// Initilize Tableau as zero
Vector T(nr * nc);
for(size_t i = 0; i < nr * nc; i++)
T[i] = 0.0;

// initialize basic variable flag as false
for(size_t j = 0; j < nc; j++)
basic[j] = false;

// For i = 0 , ... , m-1, place the Equations
// sum_j A_{i,j} * x_j + b_i <= 0 in Tableau
na = 0; // use as index of next auxillary variable
for(size_t i = 0; i < ne; i++)
{     if( b[i] > 0.0)
{     // convert to - sum_j A_{i,j} x_j - b_i >= 0
for(size_t j = 0; j < nx; j++)
T[i * nc + j] = - A[i * nx + j];
// slack variable has negative coefficient
T[i * nc + (nx + i)] = -1.0;
// auxillary variable is basic for this constraint
T[i * nc + (nx + ne + na)] = 1.0;
basic[nx + ne + na]        = true;
// right hand side
T[i * nc + (nc - 1)] = b[i];
//
++na;
}
else
{     // sum_j A_{i,j} x_j + b_i <= 0
for(size_t j = 0; j < nx; j++)
T[i * nc + j] = A[i * nx + j];
//  slack variable is also basic
T[ i * nc + (nx + i) ]  = 1.0;
basic[nx + i]           = true;
// right hand side for equations
T[ i * nc + (nc - 1) ] = - b[i];
}
}
// na is back to its original value
CPPAD_ASSERT_UNKNOWN( nc == nx + ne + na + 1 );
//
// place the equation objective equation in Tablueau
// row ne corresponds to the equation z - sum_j c_j x_j = 0
// column index for z is nx + ne + na
for(size_t j = 0; j < nx; j++)
T[ne * nc + j] = - c[j];
//
// row ne+1 corresponds to the equation w - a_0 - ... - a_{na-1} = 0
// column index for w is nx + ne + na +1
for(size_t j = 0; j < na; j++)
T[(ne + 1) * nc + (nx + ne + j)] = -1.0;
//
// fix auxillary objective so coefficients in w
// for auxillary variables are zero
for(size_t k = 0; k < na; k++)
{     size_t ja  = nx + ne + k;
size_t ia  = ne;
for(size_t i = 0; i < ne; i++)
{     if( T[i * nc + ja] != 0.0 )
{     CPPAD_ASSERT_UNKNOWN( T[i * nc + ja] == 1.0 );
CPPAD_ASSERT_UNKNOWN( T[(ne + 1) * nc + ja] == -1.0 )
ia = i;
}
}
for(size_t j = 0; j < nc; j++)
T[(ne + 1) * nc + j] += T[ia * nc + j];
// The result in column ja is zero, avoid roundoff
T[(ne + 1) * nc + ja] = 0.0;
}
//
// index of current objective
size_t iobj = ne;  // original objective z
if( na > 0 )
iobj = ne + 1; // auxillary objective w
//
// simplex interations
for(size_t itr = 0; itr < maxitr; itr++)
{     // current value for xout
for(size_t j = 0; j < nx; j++)
{     xout[j] = 0.0;
if( basic[j] )
{     // determine which row of column j is non-zero
xout[j] = std::numeric_limits<double>::quiet_NaN();
for(size_t i = 0; i < ne; i++)
{     double T_ij = T[i * nc + j];
CPPAD_ASSERT_UNKNOWN( T_ij == 0.0 || T_ij == 1.0 );
if( T_ij == 1.0 )
{     // corresponding value in right hand side
xout[j] = T[ i * nc + (nc-1) ];
}
}
}
}
if( level > 1 )
if( level > 0 )
std::cout << "itr = " << itr;
if( iobj > ne )
std::cout << ", auxillary objective w = ";
else
std::cout << ", objective z = ";
std::cout << T[iobj * nc + (nc - 1)] << "\n";
}
//
// number of variables depends on objective
size_t nv = nx + ne;   // (x, s)
if( iobj == ne + 1 )
{     // check if we have solved the auxillary problem
bool done = true;
for(size_t k = 0; k < na; k++)
if( basic[nx + ne + k] )
done = false;
if( done )
{     // switch to optimizing the original objective
iobj = ne;
}
else
nv = nx + ne + na; // (x, s, a)
}
//
// determine variable with maximuim coefficient in objective row
double cmax = 0.0;
size_t jmax = nv;
for(size_t j = 0; j < nv; j++)
{     if( T[iobj * nc + j] > cmax )
cmax = T[ iobj * nc + j];
jmax = j;
}
}
// check for solution
if( jmax == nv )
{     if( iobj == ne )
{     if( level > 0 )
std::cout << "end simplex_method\n";
return true;
}
if( level > 0 )
std::cout << "end_simples_method: no feasible solution\n";
return false;
}
//
// We will increase the j-th variable.
// Determine which row will be the pivot row.
double rmin = std::numeric_limits<double>::infinity();
size_t imin = ne;
for(size_t i = 0; i < ne; i++)
{     if( T[i * nc + jmax] > 0.0 )
{     double r =     T[i * nc + (nc-1) ] / T[i * nc + jmax];
if( r < rmin )
{     rmin = r;
imin = i;
}
}
}
if( imin == ne )
{     // not auxillary objective
if( level > 0 ) std::cout
<< "end simplex_method: objective is unbounded below\n";
return false;
}
double pivot = T[imin * nc + jmax];
//
// Which variable is changing from basic to non-basic.
// Initilaize as not yet determined.
size_t basic2not = nc;
//
// Divide row imin by pivot element
for(size_t j = 0; j < nc; j++)
{     if( basic[j] && T[imin * nc + j] == 1.0 )
{     CPPAD_ASSERT_UNKNOWN( basic2not == nc );
basic2not = j;
}
T[imin * nc + j] /= pivot;
}
// The result in column jmax is one, avoid roundoff
T[imin * nc + jmax ] = 1.0;
//
// Check that we found the variable going from basic to non-basic
CPPAD_ASSERT_UNKNOWN( basic2not < nv && basic2not != jmax );
//
// convert variable for column jmax to basic
// and for column basic2not to non-basic
for(size_t i = 0; i < nr; i++) if( i != imin )
{     double r =     T[i * nc + jmax ] / T[imin * nc + jmax];
// row_i = row_i - r * row_imin
for(size_t j = 0; j < nc; j++)
T[i * nc + j] -= r * T[imin * nc + j];
// The result in column jmax is zero, avoid roundoff
T[i * nc + jmax] = 0.0;
}
// update flag for basic variables
basic[ basic2not ] = false;
basic[ jmax ]      = true;
}
if( level > 0 ) std::cout
<< "end simplex_method: maximum # iterations without solution\n";
return false;
}
} // END_CPPAD_NAMESPACE

Input File: example/abs_normal/simplex_method.omh