********************************************************************************
* lp.ado - Linear Programming (LP) and Mixed Integer Linear Programming (MILP)
* Version: 1.0.1
* Last Updated: 25 DEC 2024
* Description:
* - implements a comprehensive framework for solving LP and MILP problems.
* - supports simplex optimization, integer constraints, and branching.
* Authors: Choonjoo Lee(Corresponding), Kyongrok Lee, Byung-ihn Lim
********************************************************************************
*! version 1.2.2 25DEC2024
*! Authors: Choonjoo Lee(Corresponding), Yongbae Ji, Kyongrok Lee, Brian Lim
*! Program: Linear Programming, Mixed Integer Linear Programming
*! Purpose: Perform LP/MILP Analysis
********************************************************************************
* HISTORY:
********************************************************************************
* 2012-10-30(THU): Basic LP/MILP Model.
* 2013-01.04(FRI): Publicly released at
* https://sourceforge.net/p/deas/code/HEAD/tree/trunk/
* 2013-08-01(THU): Presented at the 2013 Stata Conference.
* 2024-12-28(SAT): Enhanced comments and structure for clarity.
********************************************************************************
********************************************************************************
* Program Dependency Graph
********************************************************************************
* The lp.ado file contains interdependent programs and Mata functions structured
* as follows:
*
* lp (Main Program)
* ├── lpmain (LP Solver)
* │ ├── mktableau (Tableau Constructor)
* │ └── _lp_phase (Mata Phase Manager)
* │ └── lp_phase (Mata Simplex Algorithm)
* └── milp (MILP Solver)
* ├── lpmain (Relaxed LP Solver)
* ├── mktableau (Add Constraints)
* └── Recursive call to milp (Branch-and-Bound)
*
********************************************************************************
********************************************************************************
* Graph Interpretation
********************************************************************************
* 1. lp Program (Main Entry):
* - Distinguishes between LP and MILP problems.
* - Calls `lpmain` for standard LP problems.
* - Calls `milp` for MILP problems with integer constraints.
*
* 2. lpmain Program:
* - Constructs the tableau using `mktableau`.
* - Solves LP problems using `_lp_phase` Mata function for Phase I and II.
*
* 3. milp Program:
* - Handles MILP using branch-and-bound technique.
* - Uses `lpmain` for relaxed LP problems.
* - Recursively calls itself for branching on integer variables.
*
* 4. mktableau Program:
* - Prepares the tableau matrix, adding slack & artificial variables as needed.
*
* 5. Mata Functions:
* - `_lp_phase`: Manages simplex algorithm phases.
* - `lp_phase`: Implements simplex iterations for optimization.
********************************************************************************
********************************************************************************
* Program Definition Order
********************************************************************************
* 1. mktableau: Constructs the tableau matrix.
* 2. lpmain: Solves LP problems in two phases (Phase I and II).
* 3. milp: Solves MILP problems using recursive branch-and-bound.
* 4. lp: Main entry point for LP and MILP problems.
* 5. Mata Functions:
* - _lp_phase: Phase Manager.
* - lp_phase: Simplex algorithm implementation.
********************************************************************************
*! version 1.0.0 30OCT2012
*! version 1.0.1 25DEC2024
// To improve code readability, addoed comments. setting custom row names.
********************************************************************************
* Make Tableau Matrix
********************************************************************************
// -----------------------------------------------------------------------------
// This program, mktableau, constructs a Tableau Matrix for Linear Programming
// (LP) problems. The tableau is a key structure in the Simplex algorithm,
// representing the optimization problem in matrix form,
// including constraints, slack variables, and artificial variables.
// -----------------------------------------------------------------------------
// Program Definition and Syntax Parsing
// Program Definition: mktableau is defined as an rclass program to
// return matrices and scalars as results.
// Syntax Parsing:
// varlist: The list of numeric variables to include in the tableau.
// opt: Optimization goal (e.g., min or max).
// rel: Relational operators (e.g., <=, >=, =) for the constraints.
// tableau: The name of the tableau matrix to construct.
program define mktableau, rclass
syntax varlist(numeric) [if] [in], opt(string) rel(varname) tableau(name)
// Matrix Initialization:
// mkmat converts the variables (varlist) and constraints into a matrix
// (tableau). Associates row names (rel) with the matrix, representing
// the relational operators (<=, >=, =).
mkmat `varlist' `if' `in', matrix(`tableau') rownames(`rel')
// Slack and Artificial Variable Management:
// Temporary Names:
// creates placeholders (r_vec, s_mat, a_mat) for row vectors,
// slack variable matrices, and artificial variable matrices.
// Relational Operators:
// iterates over the rows of the tableau, checking each constraint's
// relational operator.
tempname r_vec s_mat a_mat
local s_names = ""
local a_names = ""
local rel_values : rownames `tableau'
forvalues i = 2/`=rowsof(`tableau')' {
matrix `r_vec' = J(rowsof(`tableau'), 1, 0)
local rel_value = word("`rel_values'", `i')
// Slack Variables for <= Constraints:
// Adds a slack variable to convert the inequality into an equation.
// Updates the slack matrix (s_mat) and names.
if ("`rel_value'" == "<" || "`rel_value'" == "<=" ) {
matrix `r_vec'[`i', 1] = 1
matrix `s_mat' = nullmat(`s_mat'), `r_vec'
local s_names = "`s_names' s`=colsof(`s_mat')'"
}
// Slack and Artificial Variables for >= Constraints:
// Adds a slack variable with a negative coefficient and
// an artificial variable.
else if ("`rel_value'" == ">" || "`rel_value'" == ">=" ) {
// slcak
matrix `r_vec'[`i', 1] = -1
matrix `s_mat' = nullmat(`s_mat'), `r_vec'
local s_names = "`s_names' s`=colsof(`s_mat')'"
// artificial
matrix `r_vec'[1, 1] = 1 // coefficients of aritificial
matrix `r_vec'[`i', 1] = 1
matrix `a_mat' = nullmat(`a_mat'), `r_vec'
local a_names = "`a_names' a`=colsof(`a_mat')'"
}
// Artificial Variables for = Constraints: adds only artificial
// variables to handle equality constraints.
else if ("`rel_value'" == "=") {
// artificial
matrix `r_vec'[1, 1] = 1 // coefficients of aritificial
matrix `r_vec'[`i', 1] = 1
matrix `a_mat' = nullmat(`a_mat'), `r_vec'
local a_names = "`a_names' a`=colsof(`a_mat')'"
}
else {
di as err "not allowed value of relational. :[`rel_value'] "
exit 198 // TODO error code confirm
}
} // end of forvalues statements
// Constructing the Final Tableau
tempname ret_tableau
// #01. Initialize objective function and variables
matrix `r_vec' = J(rowsof(`tableau'), 1, 0)
matrix `r_vec'[1,1] = 1
matrix colnames `r_vec' = "z" // Objective name
matrix `ret_tableau' = `r_vec', `tableau'[1...,1..(colsof(`tableau')-1)]
// Set custom row names. If Row names not explicitly set,
// defaults set to r1, r2, etc.
local rownames ""
forvalues i = 1/`=rowsof(`ret_tableau')' {
local rownames = "`rownames' dmu`i'"
}
matrix rownames `ret_tableau' = `rownames'
// #02. Append Slack Variables: Adds slack variables if they exist
// and counts them.
if ("`s_names'" != "") {
matrix colnames `s_mat' = `s_names'
matrix `ret_tableau' = `ret_tableau', `s_mat'
return local nslacks = colsof(`s_mat') // number of slacks
}
else return local nslacks = 0
// #03. Append Artificial Variables: Adds artificial variables if
// they exist and counts them.
if ("`a_names'" != "") {
matrix colnames `a_mat' = `a_names'
matrix `ret_tableau' = `ret_tableau', `a_mat'
return local nartificials = colsof(`a_mat') // number of artificials
}
else return local nartificials = 0
// #04. Append RHS (Right-Hand Side): Appends the RHS of the constraints
// to the final tableau.
matrix `ret_tableau' = `ret_tableau', `tableau'[1...,colsof(`tableau')]
// #05. Returning the Tableau: Returns the completed tableau matrix for
// further processing.
matrix `tableau' = `ret_tableau'
end
********************************************************************************
* LP Main - Linear Programming Main
// The provided lpmain program is a Stata script that defines the primary logic
// for solving Linear Programming (LP) problems using Phase I and Phase II of
// the two-phased revised simplex method.
********************************************************************************
program define lpmain, rclass
#del ;
// Syntax Parsing defines the required and optional arguments for the
// lpmain program.
// Required arguments:
// varlist: Variables involved in the LP problem.
// rel: Relation type (<=, >=, or =) for constraints.
// rhs: Right-hand side values of the constraints.
// opt: Optimization goal (min or max).
// Optional arguments:
// tol1: Tolerance for entering or leaving basis values (default 1e-14).
// tol2: Tolerance for B-inverse calculations (default 1e-8).
// trace: Debugging flag for detailed output during execution.
syntax varlist, rel(varname) rhs(varname) opt(string)
[
tol1(real 1e-14) tol2(real 1e-8) trace
];
#del cr
// Temporary Matrix Name Creation creates a placeholder for the Tableau
// matrix used in the simplex algorithm.
tempname tableau
// Tableau Matrix Construction:
// mktableau generates the initial tableau matrix for the LP problem.
// adds slack and artificial vars as necessary based on the rel argument.
// Metadata:
// nvars: Number of decision variables in the problem.
// nslacks: Number of slack variables added for inequalities.
// nartificials: Number of artificial variables added for equalities or
// infeasible cases.
mktableau `varlist' `rhs', opt(`opt') rel(`rel') tableau(`tableau')
local nvars : list sizeof varlist // number of variables
local nslacks = r(nslacks) // number of slacks
local nartificials = r(nartificials) // number of artificials
// LP Phase I and II Execution : This section performs calculations and
// receives results using the _lp_phase Mata code included
// in the ldea.mlib library. The following "Start of the MATA Definition
// Area" code section establishes the connection to MATA and defines
// its functionality.
// Runs the simplex algorithm in Mata for Phase I (removal of artificial
// variables) and Phase II (optimization of the objective function).
// Inputs:
// tableau: The LP tableau matrix. opt: Optimization goal (min or max).
// nvars, nslacks, nartificials: Metadata for problem structure.
// tol1, tol2: Precision tolerances. trace: Debugging output flag.
mata: _lp_phase("`tableau'", "`opt'", ///
`nvars', `nslacks', `nartificials', ///
`tol1', `tol2', "`trace'")
// Result Return:
// Returns results for use by other programs or the user.
// Returned Values:
// nvars: Number of decision variables.
// nslacks: Number of slack variables.
// nartificials: Number of artificial variables.
// tableau: The final tableau matrix after simplex iterations.
// lprslt: The final solution, including optimal values of decision
// variables and the objective function.
return local nvars = `nvars'
return local nslacks = `nslacks'
return local narticials = `nartificials'
return matrix tableau = `tableau'
return add // The return add command performs the actual task of
// adding r(lprslt) to the return list.r(lprslt)
end
program define lpmain_1, rclass
#del ;
syntax varlist, rel(varname) rhs(varname) opt(string) lprslt(name)
tableau(name) vars(real) slacks(real) artificials(real)
[
intvars(varlist) tol1(real 1e-14) tol2(real 1e-8) trace
];
#del cr
mata: _lp_phase("`tableau'", "`opt'", ///
`vars', `slacks', `artificials', ///
`tol1', `tol2', "`trace'")
tempname c_lprslt // current lprslt
matrix `c_lprslt' = r(lprslt)
matrix colnames `c_lprslt' = `: colnames(`lprslt')'
matrix rownames `c_lprslt' = `: rownames(`lprslt')'
// FIXME
// di as result _n "lprslt:"
// matrix list `lprslt', noblank nohalf noheader f(%9.6g)
// di as result _n "c_lprslt:"
// matrix list `c_lprslt', noblank nohalf noheader f(%9.6g)
if ("`intvars'" != "" && `c_lprslt'[1,1] < .) { // if MILP then,
local max_varname = ""
local max_mantissa = 0
foreach varname of varlist `intvars' {
local varvalue = ///
round(`c_lprslt'[1, colnumb(`c_lprslt',"`varname'")], `tol1')
local varvalue = `varvalue' - floor(`varvalue')
if (`varvalue' > `max_mantissa') {
local max_mantissa = `varvalue'
local max_varname = "`varname'"
}
}
if ("`max_varname'" != "") { // variables is not at all integer
tempname t_tableau t_obj t_vars t_slacks t_artificials t_rhs t_st
tempname r1_lprslt r2_lprslt temp_t
local varvalue = `c_lprslt'[1, colnumb(`c_lprslt',"`max_varname'")]
preserve
qui {
set obs `=c(N)+1'
replace `max_varname' = 1 in `c(N)'
replace `rel' = ">=" in `c(N)'
replace `rhs' = ceil(`varvalue') in `c(N)'
foreach varname of varlist `varlist' {
if ("`max_varname'" != "`varname'") {
replace `varname' = 0 in `c(N)'
}
}
}
// make tableau
mktableau `varlist' `rhs', opt(`opt') rel(`rel') tableau(`t_tableau')
local r1_vars = `vars'
local r1_slacks = r(nslacks)
local r1_artificials = r(nartificials)
// make lprslt and setup lprslt colnames and rownames
matrix `r1_lprslt' = J(1, `=(1 + `vars' + `r1_slacks')', .)
matrix `temp_t' = `t_tableau'[1...,1..`=colsof(`r1_lprslt')']
matrix colnames `r1_lprslt' = `: colnames `temp_t''
matrix rownames `r1_lprslt' = "opt_val"
// call the lp main function
lpmain `varlist', rel(`rel') rhs(`rhs') opt(`opt') ///
lprslt(`r1_lprslt') tableau(`t_tableau') ///
vars(`vars') slacks(`r1_slacks') artificials(`r1_artificials') ///
intvars(`intvars') tol1(`tol1') tol2(`tol2') `trace'
// setup result of lprslt
matrix `r1_lprslt' = r(lprslt)
/*
if (`r1_lprslt'[1,1] >= .) {
break
}
*/ restore, preserve
}
else { // select lprslt because all variables are integer
if (`lprslt'[1,1] >= .) {
matrix `lprslt' = `c_lprslt'
}
else if ("`opt'" == "max") {
if (`c_lprslt'[1,1] > `lprslt'[1,1]) {
matrix `lprslt' = `c_lprslt'
}
}
else { // else if ("`opt'" == "min") {
if (`c_lprslt'[1,1] < `lprslt'[1,1]) {
matrix `lprslt' = `c_lprslt'
}
}
}
}
else if (`c_lprslt'[1,1] < .) {
matrix `lprslt' = `c_lprslt'
}
// FIXME
di as result _n "final lprslt:"
matrix list `lprslt', noblank nohalf noheader f(%9.6g)
return matrix lprslt = `lprslt'
end
********************************************************************************
* MILP - Mixed Integer Linear Programming
// This milp program defines a recursive Mixed Integer Linear Programming (MILP)
// solution process using Stata. It combines integer constraint checking with
// recursive branch-and-bound methodology to solve optimization problems.
********************************************************************************
program define milp, rclass
#del ;
// Syntax Declaration and Initialization
// Syntax defines required inputs:
// varlist: Variables involved in the optimization.
// rel: Relation between variables (e.g., <=, >=).
// rhs: Right-hand side of the constraints.
// opt: Optimization goal (min or max).
// intvars: List of integer-constrained variables.
// Optional inputs:
// cnt: Recursive depth counter (default: 0).
// tol1, tol2: Tolerances for numeric precision.
// trace: Debugging trace flag.
// Temporary Names: tempname creates placeholders for matrices and
// variables (tableau, lprslt, etc.).
syntax varlist, rel(varname) rhs(varname) opt(string) intvars(varlist)
[
cnt(integer 0) tol1(real 1e-14) tol2(real 1e-8) trace
];
#del cr
tempname tableau lprslt baseval
// #L0. Solve Initial LP Problem: Solves the relaxed LP problem without
// considering integer constraints.
// lpmain: Solves the LP problem and stores results in:
// tableau (Problem tableau matrix), lprslt(Solution matrix).
lpmain `varlist', rel(`rel') rhs(`rhs') opt(`opt') ///
tol1(`tol1') tol2(`tol2') `trace'
matrix `tableau' = r(tableau)
matrix `lprslt' = r(lprslt)
// Debugging Output: Displays input values and results for the current iteration
// (L indicates recursion depth).
di as result _n(2) "MILP L`cnt' Input Values:"
list
matrix list `tableau', noblank nohalf noheader f(%9.6g)
di as result _n(2) "MILP L`cnt' Results: options(`opt')"
matrix list `lprslt', noblank nohalf noheader f(%9.6g)
di as text _n "--------------------------------------------------"
di as text _n
// Feasibility Check: If the LP problem is infeasible, it terminates and
// returns the results of lpmain.
if (`lprslt'[1,1] >= .) {
return add // all results of lpmain
}
else {
// Integer Constraint Validation: check that all variables is an integer
// Iterates through integer variables (intvars) and checks
// if any variable is not an integer.
// If all variables satisfy integer constraints (max_mantissa is zero),
// the solution is feasible, and the program terminates.
local max_varname = ""
local max_mantissa = 0
foreach varname of varlist `intvars' {
// because tableau and lprslt are same order
local varvalue = ///
round(`lprslt'[1, colnumb(`tableau',"`varname'")], `tol1')
local mantissa = `varvalue' - floor(`varvalue')
if (`mantissa' > `max_mantissa') {
local max_mantissa = `mantissa'
local max_varname = "`varname'"
local `baseval' = `varvalue'
}
}
// if all variables is an integer
if ("`max_varname'" == "") {
return add // all results of lpmain
}
// some variables is not an integer
else {
// #L1 : Branch 1 (>= ceil(baseval)): Adds a new constraint
// (>= ceil(baseval)) to enforce an integer solution.
// Calls milp recursively to explore this branch.
preserve
qui {
set obs `=c(N)+1'
replace `max_varname' = 1 in `c(N)'
replace `rel' = ">=" in `c(N)'
replace `rhs' = ceil(``baseval'') in `c(N)'
foreach varname of varlist `varlist' {
if ("`max_varname'" != "`varname'") {
replace `varname' = 0 in `c(N)'
}
}
}
// recursive call
milp `varlist', rel(`rel') rhs(`rhs') opt(`opt') cnt(`=`cnt'+1') ///
intvars(`intvars') tol1(`tol1') tol2(`tol2') `trace'
matrix `tableau' = r(tableau)
matrix `lprslt' = r(lprslt)
local nvars = r(nvars)
local nslacks = r(nslacks)
local nartificials = r(nartificials)
// #L2: Branch 2 (<= floor(baseval)): Adds a new constraint
// <= floor(baseval)) and recursively explores this branch.
restore, preserve
qui {
set obs `=c(N)+1'
replace `max_varname' = 1 in `c(N)'
replace `rel' = "<=" in `c(N)'
replace `rhs' = floor(``baseval'') in `c(N)'
foreach varname of varlist `varlist' {
if ("`max_varname'" != "`varname'") {
replace `varname' = 0 in `c(N)'
}
}
}
// recursive call
milp `varlist', rel(`rel') rhs(`rhs') opt(`opt') cnt(`=`cnt'+2') ///
intvars(`intvars') tol1(`tol1') tol2(`tol2') `trace'
// #L1 and #L2 are infeasible or feasible
// if #L1 is infeasible or #L2 > #L1 then select #L2
// Result Selection and Return: Compares results from Branch 1 and
// Branch 2. Retains the best feasible solution based on the
// optimization goal (min or max).
tempname L2
matrix `L2' = r(lprslt)
if ("`opt'" == "max") {
if (`lprslt'[1,1] >= . | `L2'[1,1] > `lprslt'[1,1]) {
// Select best result
matrix `tableau' = r(tableau)
matrix `lprslt' = r(lprslt)
local nvars = r(nvars)
local nslacks = r(nslacks)
local nartificials = r(nartificials)
}
}
else { // else if ("`opt'" == "min") {
if (`lprslt'[1,1] >= . | `L2'[1,1] < `lprslt'[1,1]) {
// Select best result
matrix `tableau' = r(tableau)
matrix `lprslt' = r(lprslt)
local nvars = r(nvars)
local nslacks = r(nslacks)
local nartificials = r(nartificials)
}
}
restore
// return the final results
return matrix tableau = `tableau'
return matrix lprslt = `lprslt'
return local nvars = `nvars'
return local nslacks = `nslacks'
return local narticials = `nartificials'
}
}
end
********************************************************************************
* Main Program Entry - Entry point for LP and MILP problems.
// Main program for Linear Programming and Mixed Integer Linear Programming.
// Differentiates between LP and MILP problems and delegates appropriately.
********************************************************************************
// Remove any existing program named "lp" before defining it again
// This prevents conflicts with previously defined programs
capture program drop lp
// Define a user-defined program named "lp"
// The "rclass" option allows the program to save and return results
program define lp, rclass
// Ensure compatibility with Stata version 11.0 or higher
version 11.0
// syntax checking and validation-----------------------------------------------
// rel - relational
// rhs - right hand side
// example:
// lp x1 x2 x3, min
// lp x1 x2 x3, min rel(rel_var) rhs(rhs_var)
// -----------------------------------------------------------------------------
// returns 1 if the first non-blank character in the local macro `0' is a comma,
// or if `0' is empty.
// Check if the first line of the command is empty or incorrect
// The "replay()" function validates the input
if replay() {
dis as err "vars required."
exit 198
}
// Set the command delimiter to a semicolon.
#del ;
syntax varlist(min=1) [if] [in] [using/]
[,
REL(varname) // Relational variable (default: "rel")
RHS(varname) // Right-hand side variable (default: "rhs")
MIN // Minimize optimization objective
MAX // Maximize optimization objective
INTVARS(varlist) // Integer variables (e.g., mixed-integer conditions)
ROWNAMES(varname) // specifies an optional variable in the dataset that
//contains custom row names.
TOL1(real 1e-14) // Tolerance for entering or leaving values
TOL2(real 1e-8) // Tolerance for inverse matrix calculations
TRACE // Enable/disable logging
SAVing(string) // Specify result data file name
REPLACE // Allow overwriting of result data file
];
#del cr // Restore the default delimiter to newline
// Set default value for rel to "rel" if not specified
if ("`rel'" == "") local rel = "rel"
// Set default value for rhs to "rhs" if not specified
if ("`rhs'" == "") local rhs = "rhs"
// Check if the optimization option is valid (must be either MIN or MAX)
local opt = "`min'`max'"
// Check if the optimization option is valid (must be either MIN or MAX)
if (!("`opt'" == "min" || "`opt'" == "max")) {
display as error "Error: Optimization must be MIN or MAX exclusively."
exit 198
}
// Load the dataset if the "using" option is provided
if ("`using'" != "") use "`using'", clear
// Check if a dataset is loaded and contains data
if (~(`c(N)' > 0 & `c(k)' > 0)) {
display as error "Error: A valid dataset is required!"
exit 198
}
// end of syntax checking and validation ---------------------------------------
// Disable the "more" prompt to allow smooth execution of commands
set more off
// Close any existing log file named "lp_log" if open
capture log close lp_log
// Start a new log file named "lp.log", replacing any existing file
log using "lp.log", replace text name(lp_log)
// Preserve the current dataset to allow changes w/o permanently modifying it
preserve
// Check if either an 'if' condition or an 'in' range is specified
if ("`if'" != "" | "`in'" != "") {
// Apply filtering to keep observations that satisfy the specified 'if'
// and/or 'in' conditions
qui keep `in' `if'
}
// qui (quietly) Suppresses any output or messages that would typically be
// displayed in the Stata console.
// This ensures the filtering operation is performed silently without
// interrupting the user's workflow.
// keep used to retain only the obs that satisfy certain conditions.
// Observations that do not meet the conditions are permanently removed
// from the current dataset (unless preserve was used earlier).
// in secifies a range of observations (e.g., rows or indices) to keep.
// For example, in 1/100 keeps only the first 100 observations.
// if specifies a logical condition that must be true for an observation
// to be retained. For example, if age > 18 keeps only those observations
// where the age variable is greater than 18. Combination(in and if) both
// in and if can be used together.
// -------------------------------------------------------------------------
// LP Start :
// Check if the problem is LP or MILP based on the presence of intvars.
// Solve the appropriate problem using lpmain or milp.
// Retrieve results (tableau, lprslt) and
// store metadata (nvars, nslacks, nartificials).
// Format lprslt with appropriate column and row names for clear output.
// -------------------------------------------------------------------------
// Check if the problem is LP or MILP
if ("`intvars'" == "") {
lpmain `varlist', rel(`rel') rhs(`rhs') opt(`opt') ///
tol1(`tol1') tol2(`tol2') `trace'
} // Solve as standard LP
else {
milp `varlist', rel(`rel') rhs(`rhs') opt(`opt') ///
intvars(`intvars') tol1(`tol1') tol2(`tol2') `trace'
} // If intvars is not empty, it is treated as an MILP problem.
// The milp program is called to solve the MILP problem with
// the integer constraints provided in intvars.
// Store result matrices and metadata
// tempname creates temporary names for matrices and variables.
tempname tableau lprslt temp_t
// The result matrices from the program (lpmain or milp) are saved in
// tableau: The tableau matrix representing the LP/MILP problem.
// lprslt: The resulting matrix with the optimal solution.
matrix `tableau' = r(tableau)
matrix `lprslt' = r(lprslt)
// Local Variables nvars, nslacks, nartificials store metadata:
// nvars-Number of decision variables, nslacks-Number of slack variables,
// nartificials-Number of artificial variables.
local nvars = r(nvars)
local nslacks = r(nslacks)
local nartificials = r(nartificials)
// Format results for output
// temp_t extracts a submatrix from tableau that matches the size of lprslt.
matrix `temp_t' = `tableau'[1...,1..`=colsof(`lprslt')']
// colnames copies the column names from temp_t to lprslt.
matrix colnames `lprslt' = `: colnames `temp_t''
// rownames sets the row name of lprslt to "opt_val", indicating the row
// contains the optimal value results.
matrix rownames `lprslt' = "opt_val"
// -------------------------------------------------------------------------
// REPORT:
// This section of the code is responsible for displaying the results of
// the Linear Programming (LP) or Mixed Integer Linear Programming (MILP)
// process, returning key matrices and metadata, and cleaning up the
// session by restoring the dataset and closing the log.
// -------------------------------------------------------------------------
// di as result displays the text "Input Values:" in the results window.
// _n(2) ensures a blank line above and below the message
// for better readability.
di as result _n(2) "Input Values:"
// matrix list lists the tableau matrix, which contains the input tableau
// (constraints, variables, and objective function).
// Options: noblank-Removes blank lines in the output for compactness.
// nohalf-Ensures all values are fully displayed, avoiding truncation.
// noheader: Omits the matrix name and dimensions from the display.
// f(%9.6g): Formats the matrix elements to 6 decimal places.
matrix list `tableau', noblank nohalf noheader f(%9.6g)
// di as result: Displays the text "LP Results: options(`opt')",
// indicating the optimization results (e.g., "min" or "max").
di as result _n(2) "LP Results: options(`opt')"
// matrix list: Lists the lprslt matrix, which contains the results of the
// optimization process (e.g., optimal values of variables).
// The same display options as above ensure consistent formatting.
matrix list `lprslt', noblank nohalf noheader f(%9.6g)
// di as text: adds a blank line for better visual separation in the output.
di as text _n(2) ""
// Purpose: Saves the results as retrievable return values in Stata.
// These values can be accessed using return list or programmatically by
// other commands or scripts.
// Returned Values: Matrices-tableau(Input tableau matrix), lprslt(Results
// matrix, e.g., optimal variable values).
// Scalars: nvars(Number of variables), nslacks(Number of slack variables),
// nartificials(Number of artificial variables).
return matrix tableau = `tableau'
return matrix lprslt = `lprslt'
return local nvars = `nvars'
return local nslacks = `nslacks'
return local narticials = `nartificials'
// set more on: Re-enables the "more" prompt, which was disabled earlier
// (set more off).
// restore, preserve: Restores the original dataset to its previous state,
// undoing any changes made during the program execution.
// The preserve option ensures that the restoration does not overwrite
// anything accidentally.
// log close lp_log: Closes the log file lp_log and saves it.
// This ensures that all output is properly recorded.
set more on
restore, preserve
log close lp_log
end
********************************************************************************
* Mata Definition Area
********************************************************************************
// Start of the MATA Definition Area -------------------------------------------
// The code defines two key functions in Mata for solving Linear Programming
// (LP) problems:
// _lp_phase: Manages the overall process for Phase I and Phase II of the
// Simplex algorithm. lp_phase: Performs the actual simplex iterations,
// including variable classification and optimization.
version 10
mata:
mata set matastrict on
// _lp_phase function sets up and controls the LP optimization process.
void function _lp_phase (
string scalar tableau,
string scalar opt,
real scalar vars,
real scalar slacks,
real scalar artificials,
real scalar tol1,
real scalar tol2,
string scalar trace )
{
real matrix M, VARS
real fcols
struct BoundCond matrix boundM
struct LpParam scalar param
struct LpResultStruct scalar lpresult
// 1st. Load Matrix and Variables: Loads the input tableau (M) from Stata.
// Defines variable indices for the tableau.
M = st_matrix(tableau)
VARS = (0, 1..vars+slacks, -1..-artificials, 0)
// 2rd. Set Boundaries for Variables: Defines default boundaries [0,∞].
// 0 <= weight, slacks, atrificials <= INFINITE
boundM = J(1, cols(M), BoundCond());
for (i=1; i<cols(M); i++) {
boundM[1,i].val = 0; boundM[1,i].lower = 0; boundM[1,i].upper = .
}
// 3th. Parameter Initialization: Initializes parameters based on user
// input, including Optimization goal(min or max), tolerances(tol1, tol2).
// debugging trace flag.
param.minYn = (opt == "min"); // 0: max, 1: min
param.vars = vars
param.slacks = slacks
param.artificials = artificials
param.tol1 = tol1
param.tol2 = tol2
param.trace = trace
param.tracename = "LP for RSM"
// Call lp_phase: executes the simplex algorithm, handling both Phase I
// (artificial variables) and Phase II (objective optimization).
lpresult = lp_phase(M, boundM, VARS, param)
// -------------------------------------------------------------------------
// Process Results:
// Returns: Final variable values. Objective function values. Stores
// results in the r(lprslt) matrix in Stata.
// -------------------------------------------------------------------------
if(lpresult.rc) {
LPRSLT = J(1, 1+param.vars+param.slacks, .)
}
else {
// lpresult = theta(1) + vars + slacks
LPRSLT = J(1, param.vars+param.slacks, 0)
for (j=1; j<=rows(lpresult.XB) ; j++) {
if (VARS[1,j+1] > 0) LPRSLT[1, VARS[1,j+1]] = lpresult.XB[j, 1]
}
LPRSLT = lpresult.xVal, LPRSLT
}
if (param.trace == "trace") {
msg = sprintf("%s-FINAL", param.tracename);
// printf("\n%s: original VARS.\n", msg); orgVARS
printf("\n%s: VARS.\n", msg); VARS
printf("\n%s: XB.\n", msg); lpresult.XB
printf("\n%s: LPRSLT.\n", msg); LPRSLT
}
st_matrix("r(lprslt)", LPRSLT)
}
/**
* @param VARS - Variable Index Matrix
* [z, B, N, b]'s index in the original Tableau
* @param M - Tableau: [z, A, S, Af, b] --> [z, B, N, b]
* @param phase - if have artificials, then phase 1 and 2,
* otherwise only phase 2
* @param param - parameter struct for Lp RSM
*
* @return result of LP
*/
struct LpResultStruct function lp_phase (
real matrix M,
struct BoundCond matrix boundM,
real matrix VARS,
struct LpParam scalar param )
{
real scalar phase, mrows, mcols, j, idx
string scalar tracename
real vector reorderidx, bfsidx, nonbfsidx
real vector coef_of // coefficient of objective function
struct LpParamStruct scalar lpParam
struct LpResultStruct scalar lpResult
// lp_phase function performs the core simplex algorithm steps.
// Validation Checks: ensures valid input for the optimization goal.
if (param.minYn >= .) { //
displayas("err");
_error(3351, "You have to set the minimization(1) or maximization(0) "
+ "at the LpParam.minYn")
}
// Initialize Objective Function: Extracts and temporarily removes the
// coefficients of the objective function from the tableau.
coef_of = M[1, 2..1+param.vars] // keep the objective function
replacesubmat(M, 1, 2, J(1, param.vars, 0))
// initialize matrix.
if (param.trace == "trace") {
displayas("txt")
printf("\n\n%s: initialize matrix.\n", param.tracename); M
}
mrows = rows(M); mcols = cols(M)
// Classify Variables: Identifies basic (bfsidx) and non-basic (nonbfsidx)
// variables based on tableau structure.
bfsidx = J(1, mrows-1, .); nonbfsidx = J(1, 0, .)
for (j = 2+param.vars; j <= mcols-1; j++) {
T = M[2::mrows,j]
if ((sum(T :!= 0) == 1) && (sum(T) == 1)) {
maxindex(T, 1, i, w); bfsidx[i] = j
}
else nonbfsidx = nonbfsidx, j
}
reorderidx = (1, bfsidx[1,], 2..1+param.vars, nonbfsidx[1,], mcols)
VARS = VARS[,reorderidx];
M = M[,reorderidx]; boundM = boundM[,reorderidx]
if (param.trace == "trace") {
displayas("txt")
printf("\n%s: classify basic and nonbasic.\n", tracename); M; VARS
}
// set the lp's parameters
lpParam.dmus = param.vars
lpParam.slacks = param.slacks
lpParam.artificials = param.artificials
lpParam.tol1 = param.tol1
lpParam.tol2 = param.tol2
lpParam.trace = param.trace
// Phase I - Artificial Variables: Solves the problem in Phase I to
// eliminate artificial variables.
if (param.artificials > 0) {
phase = 1
lpParam.minYn = 1; // min because of phase 1
tracename = param.tracename + "-PI"
lpResult = lp(M, boundM, VARS, 0, phase, tracename, lpParam)
if (lpResult.rc) return(lpResult)
}
// Phase II - Objective Optimization: Restores the original objective
// function and solves the optimization problem in Phase II.
phase = 2
lpParam.minYn = param.minYn // according to the optimization.
tracename = param.tracename + "-PII"
// set the objective function.
mcols = cols(M)
for (j=2; j<mcols; j++) {
idx = VARS[1,j]
if (0 < idx && idx <= param.vars) {
M[1,j] = coef_of[idx] // according to variable's index
}
}
lpResult = lp(M, boundM, VARS, 0, phase, tracename, lpParam)
// return result.
return(lpResult)
}
end
// End of the MATA Definition Area ---------------------------------------------
********************************************************************************
* Copyright and Licensing Information
********************************************************************************
/**
* Program: lp.ado (linear programming, mixed-integer programming Analysis)
* Version: 1.2.2 (28DEC2024)
* Copyright (c) 2009-2024 by Choonjoo Lee, Rep. of Korea. All Rights Reserved.
* This product includes software developed at the Lab. of Choonjoo Lee.
* Contributors: Kyung-rok Lee, Byung-ihn Lim.
*
* License:
* This program is distributed under the following license terms:
*
* Academic Free License v3.0:
* - The program may be used, modified, and redistributed freely for academic,
* educational, and research purposes.
* - Proper attribution to the authors is required in any derived work
* or publication.
*
* Commercial Restriction:
* - Commercial use, including distribution, modification, or sale for profit,
* is strictly prohibited.
* - For commercial licensing, please contact the authors directly.
*
* Warranty Disclaimer:
* - This software is provided "as-is" without warranty of any kind.
* - The authors are not responsible for any damage or data loss resulting
* from the use of this software.
*
* Citation:
* If you use this program in your research, please cite it as follows:
* Lee, C., Lee, K., and Lim, B. (2024). "LP/MILP Analysis: A Stata ADO Program.
*
* ---------------------------------------------------------------------------
* For inquiries, please contact: Choonjoo Lee <dea.stata@gmail.com>
* ---------------------------------------------------------------------------
*/