R/linp.R In limSolve: Solving Linear Inverse Models

Documented in linp

```##==============================================================================
## linp        :
## Solves a linear programming problem,
## Minimise               Cost
## subject to             E*x=f
##                        G*x>=h
##
## Note: uses lp from package lpSolve
## This R-code sometimes fails and terminates R
##      for very small problems that are repeated frequently...
##==============================================================================

linp <- function(E=NULL, F=NULL, G=NULL, H=NULL, Cost,
ispos=TRUE, int.vec=NULL, verbose=TRUE, ...)  {

## input consistency
if (! is.matrix(E) & ! is.null(E))
E <- t(as.matrix(E))
if (! is.matrix(G) & ! is.null(G))
G <- t(as.matrix(G))

## problem dimension
Neq  <- nrow(E)   # Number equalities
Nx   <- ncol(E)   # Number unknowns
Nin  <- nrow(G)   # Number inequalities

if (is.null(Nx )) Nx  <- ncol(G)
if (is.null(Nin)) Nin <- 0
if (is.null(Neq)) Neq <- 0

NumEq <- Nin+Neq  # total number of equations

## consistency of input
if (!is.null(G)) {
if (ncol(G)   != Nx)
stop("cannot solve linear programming problem - E and G not compatible")
if (length(H) != Nin)
stop("cannot solve linear programming problem - G and H not compatible")
}

if (!is.null(E)) {
if (length(F) != Neq)
stop("cannot solve linear programming problem - E and F not compatible")
}

if (length(Cost)!= Nx)
stop("cannot solve linear programming problem - Cost not compatible")

IsError <- FALSE

## con: constraints ; rhs: right hand side

## the equalities:
con   <- E
rhs   <- F
dir   <- rep("==",Neq)

## inequalities:
if (Nin > 0)  {
con   <- rbind(con,G)
rhs   <- c(rhs,H)
dir   <- c(dir,rep(">=",Nin))
}

if (!ispos)  {
con  <- cbind(con, -con)
Cost <- c(Cost, -Cost)
if (! is.null(int.vec))
int.vec<-c(int.vec,int.vec+Nx)
}

## the solution
sol    <- lp("min",Cost,con,dir,rhs,int.vec=int.vec,...)
mode   <- sol\$status
## from the original code lp_lib.h
if (mode == -5) {
print("unknown error")
IsError<-TRUE
}
else if (mode == -4) {
print("data ignored")
IsError<-TRUE
}
else if (mode == -3) {
print("no bfp")
IsError<-TRUE
}
else if (mode == -2) {
print("no memory")
IsError<-TRUE
}
else if (mode == -1) {
print("problem not run")
IsError<-TRUE
}
else if (mode == 2) {
print("problem infeasible")
IsError<-TRUE
}
else if (mode == 3) {
print("problem unbounded")
IsError<-TRUE
}
else if (mode == 4) {
print("problem degenerate")
IsError<-TRUE
}
else if (mode == 5) {
print("problem failed")
IsError<-TRUE
}
else if (mode == 6) {
print("user aborted")
IsError<-TRUE
}
else if (mode == 7) {
print("timed out")
IsError<-TRUE
}
else if (mode == 8) {
print("running")
IsError<-TRUE
}
else if (mode == 9) {
print("presolved")
IsError<-TRUE
}

X           <- sol\$solution
if (!ispos)
X <- X[1:Nx]-X[(Nx+1):(2*Nx)]
solutionNorm<- sol\$objval

## Total residual norm
residual <- 0
if (!is.null(E))
residual <- sum(abs(E %*% X - F))+residual
if (!is.null(G))   {
ineq     <- G %*% X - H
residual <- residual -sum(ineq[ineq<0])
}
xnames <- colnames(E)
if (is.null(xnames))
xnames <- colnames(G)
if (is.null(xnames))
xnames <- names(Cost[1:Nx])
names (X) <- xnames

return(list(X=X,                        # vector containing the solution of the linear programming problem.
residualNorm=residual,      # scalar, the sum of residuals of equalities and violated inequalities
solutionNorm=solutionNorm,  # scalar, the value of the minimised cost function
IsError=IsError,            # if an error occurred
type="linp"))

}
```

Try the limSolve package in your browser

Any scripts or data that you put into this service are public.

limSolve documentation built on Aug. 14, 2017, 3:01 p.m.