Nothing
solveLP <- function( cvec, bvec, Amat, maximum=FALSE,
const.dir = rep( "<=", length( bvec ) ),
maxiter=1000, zero=1e-9, tol=1e-6, dualtol = tol,
lpSolve=FALSE, solve.dual=FALSE, verbose = 0 )
{
result <- list() # list for results that will be returned
result$status <- 0
rdigits <- -round( log10( zero ) )
nVar <- length(cvec) # number of variables
nCon <- length(bvec) # number of constraints
if( !all.equal( dim( Amat ), c( nCon, nVar ) ) == TRUE ) {
stop( paste( "Matrix A must have as many rows as constraints (=elements",
"of vector b) and as many columns as variables (=elements of vector c).\n" ) )
}
if( length( const.dir ) != nCon ) {
stop( paste( "'const.dir' must have as the elements as constraints",
"(=elements of vector b).\n" ) )
}
if( sum( const.dir == ">=" | const.dir == ">" | const.dir == "=" |
const.dir == "==" | const.dir == "<=" | const.dir == "<" ) < nCon ) {
stop( "'const.dir' may only contain '>=', '>', '=', '==', '<=' or '<'" )
}
if( any( const.dir %in% c( "=", "==" ) ) && ( ! lpSolve ) ) {
warning( "solveLP() might return incorrect results",
" if the model includes equality constraints",
" and argument 'lpSolve' is 'FALSE';",
" please check if solveLP() returns the same results",
" with argument 'lpSolve' equal to 'TRUE';",
" more information on this bug available at",
" linprog's R-Forge site" )
}
## Labels
if( is.null(names(cvec))) {
clab <- as.character(1:nVar)
} else {
clab <- names(cvec)
}
if( is.null(names(bvec))) {
blab <- as.character(1:nCon)
} else {
blab <- names(bvec)
}
const.dir2 <- rep( 0, nCon )
const.dir2[ const.dir == ">=" | const.dir == ">" ] <- 1
const.dir2[ const.dir == "<=" | const.dir == "<" ] <- -1
## lpSolve
if( lpSolve ) {
if( maximum ) {
direction <- "max"
} else {
direction <- "min"
}
lpres <- lp( direction = direction, cvec, Amat, const.dir, bvec )
result$lpStatus <- lpres$status
if( result$lpStatus == 0 ) {
if( min( lpres$solution ) < -tol ) {
result$lpStatus <- 7
} else if( max( ( bvec - c( Amat %*% lpres$solution ) ) *
const.dir2 ) > tol ) {
result$lpStatus <- 3
}
}
if( result$lpStatus != 0 ) {
result$status <- 1
} else {
result$solution <- lpres$solution
names( result$solution ) <- clab
result$opt <- lpres$objval
## Results: Constraints
result$con <- data.frame( actual=NA, dir=const.dir, bvec=bvec, free=NA )
result$con$actual <- round( c( Amat %*% result$solution ), digits=rdigits )
names( result$con$actual ) <- blab
result$con$free <- round( result$con$bvec - result$con$actual, digits=rdigits )
result$con$free[ const.dir2 == 1 ] <- -result$con$free[ const.dir2 == 1 ]
result$con$free[ const.dir2 == 0 ] <- -abs( result$con$free[ const.dir2 == 0 ] )
}
} else {
## Simplex algorithm
iter1 <- 0
iter2 <- 0
## Slack Variables
for(i in 1:nCon) clab <- c( clab, paste("S", blab[i] ) )
cvec2 <- c( cvec, rep( 0, nCon ) )
## Tableau ( Basic Variables, Slack,Variables, P0, Z-C )
Tab <- rbind( cbind( -Amat * const.dir2, diag( 1, nCon, nCon ), -bvec * const.dir2 ),
c( cvec2 * (-1)^maximum, 0 ) )
rownames(Tab) <- c( blab, "Z-C" )
colnames(Tab) <- c( clab, "P0" )
if( verbose >= 3 ) {
print("initial Tableau")
print(Tab)
}
## searching for feasible solution for starting
# basis: Zero Solution ( Basic Variables = Slack Variables )
basis <- c( (nVar+1) : (nVar+nCon) )
if(min(Tab[ 1:nCon, nVar+nCon+1]) < 0 ) {
Tab2 <- Tab
Tab2 <- rbind( Tab2, rep(0, ncol(Tab2) ) )
rownames(Tab2)[nCon+2] <- "M Z-C" # additional artificial 'Z-C' row
basis2 <- basis
nArt <- 0 # number of artificial variables
for(i in 1:nCon) {
if(Tab[ i, nVar+nCon+1] < 0 ) {
Tab2[ i, ] <- -Tab2[ i, ]
Tab2 <- cbind( Tab2[ , 1:(nVar+nCon+nArt) ], rep(0,nCon+2),
Tab2[ , (nVar+nCon+nArt+1)] )
nArt <- nArt + 1
colnames(Tab2)[ nVar+nCon+nArt ] <- paste("M", rownames(Tab2)[i] )
Tab2[ i, nVar+nCon+nArt ] <- 1
Tab2[ nCon+2, nVar+nCon+nArt ] <- 1
# put artificial variables in basis
rownames(Tab2)[ i ] <- paste("M", rownames(Tab2)[i] )
basis2[i] <- nVar+nCon+nArt
}
}
for(i in 1:nCon) { # artificial objective function (Z-C)
if(Tab[ i, nVar+nCon+1] < 0 ) {
Tab2[nCon+2, 1:(nVar+nCon+nArt)] <- Tab2[nCon+2, 1:(nVar+nCon+nArt)] -
Tab2[ i , 1:(nVar+nCon+nArt)]
}
}
for(i in 1:nCon) { # value of artificial objective function
Tab2[nCon+2, nVar+nCon+nArt+1 ] <- Tab2[nCon+2, nVar+nCon+nArt+1 ] -
Tab2[i, nVar+nCon+nArt+1] * Tab2[nCon+2, basis[i] ]
}
colnames(Tab2)[ nVar+nCon+nArt+1 ] <- "P0"
if( verbose >= 3 ) {
print("initial Tableau for Phase 1")
print(Tab2)
}
## Simplex algorithm (Phase 1)
while( min( Tab2[ nCon+2, 1:(nVar+nCon+nArt) ] ) < -zero & iter1 < maxiter) {
iter1 <- iter1 + 1
## Pivot
Tab[ abs(Tab) < zero ] <- 0
# pcolumn <- which.min( Tab2[ nCon+2, 1:(nVar+nCon+nArt) ]) # Pivot column
decval <- array( NA, nVar+nCon )
for( pcolumnt in 1:(nVar+nCon+nArt) ) {
if( Tab2[ nCon+2, pcolumnt ] < 0 ) {
rwerte <- Tab2[ 1:nCon, nVar+nCon+nArt+1 ] / Tab2[ 1:nCon , pcolumnt ]
# R-values
rwerte[ Tab2[1:nCon, pcolumnt ] <= 0 ] <- max(rwerte,na.rm=TRUE)+1
prow <- which.min( rwerte ) # Pivot row
if( length( rwerte[ !is.na(rwerte) & is.finite(rwerte) ] ) >= 1 ) {
decval[ pcolumnt ] <- Tab2[ nCon+2, pcolumnt ] *
min( rwerte[ !is.na(rwerte) & is.finite(rwerte) ] )
}
}
}
if( min( decval, na.rm=TRUE ) < -zero ) {
pcolumn <- which.min( decval ) # Pivot column
} else {
pcolumn <- which.min( Tab2[ nCon+2, 1:(nVar+nCon+nArt) ]) # Pivot column
}
rwerte <- Tab2[ 1:nCon , nVar+nCon+nArt+1 ] / Tab2[ 1:nCon , pcolumn ] # R-values
rwerte[ Tab2[1:nCon, pcolumn ] <= 0 ] <- max(rwerte, na.rm=TRUE)+1
prow <- which.min( rwerte ) # Pivot row
if( verbose >=2 ) {
cat( paste( "\nPivot Column:", as.character(pcolumn),
"(",colnames(Tab2)[pcolumn],")\n" ) )
cat( paste( "Pivot Row:", as.character( prow ),
"(", rownames(Tab2)[prow], ")\n\n" ) )
}
## New Basis
basis[prow] <- pcolumn
rownames(Tab2)[prow] <- colnames(Tab2)[pcolumn]
## New Tableau
Tab2[ prow, ] <- Tab2[ prow, ] / Tab2[ prow, pcolumn ]
for( i in 1:(nCon+2) ) {
if( i != prow ) {
Tab2[ i, ] <- Tab2[ i, ] - Tab2[ prow, ] * Tab2[ i, pcolumn ]
}
}
if( verbose >= 4 ) print(Tab2)
}
if(iter1 >= maxiter ) warning("Simplex algorithm (phase 1) did not reach optimum.")
Tab <- cbind( Tab2[ 1:(nCon+1), 1:(nCon+nVar) ],
Tab2[ 1:(nCon+1), nVar+nCon+nArt+1 ] )
if( verbose >= 3 ) {
print("New starting Tableau for Phase II")
print(Tab)
}
}
## Simplex algorithm (Phase 2)
while( min( Tab[ nCon+1, 1:(nVar+nCon) ] ) < -zero & iter2 < maxiter ) {
iter2 <- iter2 + 1
## Pivot
Tab[ abs(Tab) < zero ] <- 0
# pcolumn <- which.min( Tab[ nCon+1, 1:(nVar+nCon) ]) # Pivot column
decval <- array( NA, nVar+nCon )
for( pcolumnt in 1:(nVar+nCon) ) {
if( Tab[ nCon+1, pcolumnt ] < 0 ) {
rwerte <- Tab[ 1:nCon , nVar+nCon+1 ] / Tab[ 1:nCon , pcolumnt ]
# R-values
rwerte[ Tab[1:nCon, pcolumnt ] <= 0 ] <- max(rwerte,na.rm=TRUE)+1
prow <- which.min( rwerte ) # Pivot row
if( length( rwerte[ !is.na(rwerte) & is.finite(rwerte) ] ) >= 1 ) {
decval[ pcolumnt ] <- Tab[ nCon+1, pcolumnt ] *
min( rwerte[ !is.na(rwerte) & is.finite(rwerte) ] )
}
}
}
if( min( decval, na.rm=TRUE ) < -zero ) {
pcolumn <- which.min( decval ) # Pivot column
} else {
pcolumn <- which.min( Tab[ nCon+1, 1:(nVar+nCon) ]) # Pivot column
}
rwerte <- Tab[ 1:nCon , nVar+nCon+1 ] / Tab[ 1:nCon , pcolumn ] # R-values
rwerte[ Tab[1:nCon, pcolumn ] <= 0 ] <- max(rwerte,na.rm=TRUE)+1
prow <- which.min( rwerte ) # Pivot row
if( verbose >= 2 ) {
cat( paste( "\nPivot Column:", as.character(pcolumn),
"(",colnames(Tab)[pcolumn],")\n" ) )
cat( paste( "Pivot Row:", as.character( prow ) ,
"(",rownames(Tab)[prow],")\n\n") )
}
## New Basis
basis[prow] <- pcolumn
rownames(Tab)[prow] <- colnames(Tab)[pcolumn]
## New Tableau
Tab[ prow, ] <- Tab[ prow, ] / Tab[ prow, pcolumn ]
for( i in 1:(nCon+1) ) {
if( i != prow ) {
Tab[ i, ] <- Tab[ i, ] - Tab[ prow, ] * Tab[ i, pcolumn ]
}
}
if( verbose >= 4 ) print(Tab)
}
if(iter2 >= maxiter ) warning("Simplex algorithm (phase 2) did not reach optimum.")
## Results: Basic Variables
basvar <- matrix( NA, nCon, 1 )
colnames(basvar) <- c("opt")
rownames(basvar) <- rep("a",nCon)
for( i in 1:nCon ) {
rownames(basvar)[i] <- clab[sort(basis)[i]]
basvar[i,1] <- Tab[ which(basis==sort(basis)[i]), nVar+nCon+1 ]
}
## Results: All Variables (Including Slack Variables)
allvar <- data.frame( opt=rep( NA, nVar+nCon ), cvec=cvec2, min.c=NA,
max.c=NA, marg=NA, marg.reg=NA )
rownames(allvar) <- clab
for( i in 1:(nVar+nCon) ) {
if(i %in% basis ) {
allvar$opt[ i ] <- Tab[ which(basis==i), nVar+nCon+1 ]
## Stability of Basic Variables
quot <- Tab[ nCon+1, 1:(nVar+nCon) ] / Tab[ which(basis==i), 1:(nVar+nCon) ]
if( maximum ) {
if(max(quot[!is.na(quot)]) > 0 ) {
suppressWarnings(
allvar$min.c[ i ] <- cvec2[ i ] - min(quot[quot>0 & !is.na(quot)])
)
}
if(min(quot[!is.na(quot) & is.finite(quot)]) < 0 ) {
if(max(quot[quot<0 & !is.na(quot)]) > -1e14 ) {
allvar$max.c[ i ] <- cvec2[ i ] - max(quot[quot<0 & !is.na(quot)])
} else {
allvar$max.c[ i ] <- Inf
}
} else {
allvar$max.c[ i ] <- Inf
}
} else {
if(max(quot[!is.na(quot)]) > 0 ) {
suppressWarnings(
allvar$max.c[ i ] <- cvec2[ i ] + min(quot[quot>0 & !is.na(quot)])
)
}
if(min(quot[!is.na(quot)]) < 0 ) {
if(max(quot[quot<0 & !is.na(quot)]) > -1e14 ){
allvar$min.c[ i ] <- -cvec2[ i ] + max(quot[quot<0 & !is.na(quot)])
} else {
allvar$min.c[ i ] <- NA
}
} else {
allvar$min.c[ i ] <- NA
}
}
} else {
allvar$opt[ i ] <- 0
if( i <= nVar ) {
if( maximum ) {
allvar$min.c[ i ] <- -Inf
allvar$max.c[ i ] <- Tab[ nCon+1, i ] + cvec2[i]
} else {
allvar$min.c[ i ] <- 99#-Tab[ nCon+1, i ] - cvec2[i]
allvar$max.c[ i ] <- 77#Inf
}
}
}
allvar$cvec[ i ] <- cvec2[ i ]
# marginal contribution to objective function (Shadow prices)
if( !( ( i %in% basis ) & ( i <= nVar ) ) ) {
allvar$marg[ i ] <- Tab[ nCon+1, i ] * (-1)^maximum
if( !( i %in% basis ) & ( i > nVar ) ) {
if( maximum ) {
allvar$max.c[ i ] <- Tab[ nCon+1, i ] #* (-1)^maximum
allvar$min.c[ i ] <- -Inf
} else {
allvar$min.c[ i ] <- -Tab[ nCon+1, i ]
allvar$max.c[ i ] <- Inf
}
}
quot <- Tab[ 1:nCon , nVar+nCon+1 ] / Tab[ 1:nCon, i ]
suppressWarnings(
if( !( i %in% basis) ) {
allvar$marg.reg[ i ] <- min(quot[quot>0 & !is.na(quot)])
} else {
allvar$marg.reg[ i ] <- NA
}
)
}
}
allvar$min.c[ allvar$min.c > 1e16 ] <- Inf
allvar$min.c[ allvar$min.c < -1e16 ] <- -Inf
## Results: Constraints
con <- data.frame( actual=NA, dir=const.dir, bvec=bvec, free=NA, dual=NA, dual.reg=NA )
names( con$actual ) <- blab
for(i in 1: nCon) {
if( (i+nVar) %in% basis ) {
con$actual[ i ] <- round( bvec[i] + Tab[ which((i+nVar)==basis),
nVar+nCon+1 ] * const.dir2[ i ], digits=rdigits )
} else {
con$actual[ i ] <- round( bvec[i], digits=rdigits )
}
if( -allvar$opt[ i+nVar ] == 0 ) {
con$dual[ i ] <- allvar$marg[ i+nVar ] * (-1)^maximum
con$dual.reg[ i ] <- allvar$marg.reg[ i+nVar ]
} else {
con$dual[ i ] <- 0
con$dual.reg[ i ] <- allvar$opt[ i+nVar ]
}
}
con$free <- round( con$bvec - con$actual, digits=rdigits )
con$free[ const.dir2 == 1 ] <- -con$free[ const.dir2 == 1 ]
con$free[ const.dir2 == 0 ] <- -abs( con$free[ const.dir2 == 0 ] )
result$opt <- round( -Tab[ nCon+1, nCon+nVar+1 ], digits=rdigits ) * (-1)^maximum
result$iter1 <- iter1
result$iter2 <- iter2
result$allvar <- round( allvar, digits=rdigits )
result$basvar <- round( basvar, digits=rdigits )
result$solution <- result$allvar$opt[ 1 : nVar ]
names( result$solution ) <- clab[ 1: nVar ]
result$con <- con
if( verbose >= 1 ) result$Tab <- Tab
if( iter1 >= maxiter ) result$status <- 4
if( iter2 >= maxiter ) result$status <- 5
}
if( result$status == 0 ) {
if( min ( result$con$free ) < - tol ) {
result$status <- 3
}
}
## solving the dual problem
if( solve.dual && result$status == 0 ) {
if( any( const.dir2 == 0 ) ) {
stop( paste( "At the moment the dual problem can not be solved",
"with equality constraints" ) )
}
if( maximum ) {
const.dir.dual <- rep(">=",nVar)
} else {
const.dir.dual <- rep("<=",nVar)
}
result$con$dual.p <- result$con$dual
dualres <- solveLP( cvec = bvec * const.dir2 * (-1)^maximum, bvec = cvec,
Amat = t( Amat * const.dir2 ) * (-1)^maximum, maximum = !maximum,
const.dir = const.dir.dual, maxiter = maxiter, zero = zero,
tol = dualtol, lpSolve = lpSolve, verbose = verbose )
result$dualStatus <- dualres$status
if( result$dualStatus == 0 ) {
result$con$dual <- dualres$solution
} else {
result$status <- 2
}
}
## List of Results
result$maximum <- maximum
result$lpSolve <- lpSolve
result$solve.dual <- solve.dual
result$maxiter <- maxiter
class(result) <- "solveLP"
result
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.