R/linprog.R

In pracma: Practical Numerical Math Functions

##
##  l i n p r o g . R  Linear Programming Solver
##


linprog <- function(cc, A = NULL, b = NULL, Aeq = NULL, beq = NULL,
                    lb = NULL, ub = NULL, x0 = NULL, I0 = NULL,
                    bigM = 100, maxiter = 20, maximize = FALSE)
{
    .lp.check(cc, A, b, Aeq, beq, lb, ub, x0, I0)
    if (maximize) cc <- -cc

    # there are only inequality constraints
    if (is.null(Aeq) && all(b >= 0)) {
        # go and find a base vector and index set
        Pin <- .lp.create(cc, A, b, lb, ub)
        Sin <- .lp.solve(Pin$cc, Pin$A, Pin$b, Pin$x0, Pin$I0,
                         maxiter = maxiter)

    # find a feasible base vector first
    } else {
        # there are only equality constraints
        if (is.null(A) && is.null(lb) && is.null(ub)) {
            A <- Aeq
            b <- beq
            cc0 <- cc

        # there are both equality and inequality constraints
        } else {
            Pin <- .lp.create(cc, A, b, lb, ub)
            Ain <- Pin$A
            bin <- Pin$b
            cc0 <- Pin$cc
            if (!is.null(Aeq)) {
                Aeq <- cbind(Aeq, matrix(0, nrow = nrow(Aeq),
                                    ncol = ncol(Ain)-ncol(Aeq)))
                A <- rbind(Ain, Aeq)
                b <- c(Pin$b, beq)
            } else {
                A <- Ain
                b <- bin
            }
        }

        # only equalities; now guarantee that b >= 0
        inds <- which(b < 0)
        if (length(inds) > 0) {
            A[inds, ] <- -A[inds, ]
            b[inds] <- -b[inds]

        }
        # all constraints combined in one equality matrix A and vectors b, cc
        m <- nrow(A); n <- ncol(A)
        A0 <- cbind(A, diag(m))
        b0 <- b

        x0 <- c(rep(0, n), b0)
        I0 <- (n+1):(n+m)

        if (bigM > 0) {                 # big-M approach
            M  <- bigM
            k <- 1; kmax <- 8
            while (k < kmax) {
                f <- c(cc0, rep(M, m))
                Sin <- .lp.solve(f, A0, b0, x0, I0, maxiter = maxiter)
                if (Sin$errno < 0 || all(Sin$x[I0] == 0)) break
                M <- 10 * M
                k <- k + 1
            }
            # browser()
            if (k >= kmax) {
                Sin$errno <- -4
            }
        } else {                        # phase-I approach
            f <- rep(1, n+m)
            Sin <- .lp.solve(f, A0, b0, x0, I0, maxiter = maxiter)
            if (Sin$errno < 0) {
                Sin$errno <- -5
            } else {
                x1 <- Sin$x[1:n]; z1 <- Sin$x[(n+1):(n+m)]
                I1 <- which(x1 != 0)
                if (length(I1) != m) {
                    Sin$errno <- -5
                } else if (all(z1 == 0)) {
                    Sin <- .lp.solve(cc0, A, b, x1, I1, maxiter = maxiter)
                } else
                    Sin$errno <- -2
            }
        }
    }
    errno <- Sin$errno
    errmsg <- c("Maximum no. of iterations exceeded.",
                "Problem is most likely unfeasible.",
                "Problem is most likely unbounded.",
                "Big-M appraoch not successful.",
                "Phase-I approach not successful.")

    if (errno < 0) {
        mess <- errmsg[abs(errno)]
        return(list(x = NA, fval = NA, errno = errno, message = mess))
    } else {
        mess <- "Solver LP converged successfully."
        x <- Sin$x[1:length(cc)]
        y <- Sin$y[1:length(b)]
        f <- sum(cc * x)
        if (maximize) f <- -f
        return(list(x = x,  fval = f, errno = errno, message = mess))
    }
    cat("LP solver stopped.\n")
}


.lp.check <- function(cc, A, b, Aeq, beq, lb, ub, x0, I0) {
    stopifnot(is.numeric(cc))
    if (is.numeric(lb) && length(lb) != length(cc) ||
        is.numeric(ub) && length(ub) != length(cc))
        stop("Arguments 'cc', 'lb', and 'ub' must have equal lengths.")

    if (!is.null(A)) {
        if (is.null(b)) stop("Vector 'b' must be present if 'A' is.")
        if (!is.matrix(A)) stop("Argument 'A' must be a matrix.")
        m <- nrow(A); n <- ncol(A)
        if (length(cc) != n) stop("Length of 'cc' does not fit to 'ncol(A)'.")
        if (length(b) != m)  stop("Length of 'b' does not fit to 'nrow(A)'.")
    }
}


.lp.create <- function(cc, Ain, bin, lb, ub) {
    A <- Ain
    b <- bin
    if (is.null(Ain)) {
        m_in <- 0
        n <- length(cc)
    } else {
        m_in <- nrow(Ain)
        n <- ncol(Ain)
    }

    # upper bounds
    if (!is.null(ub)) {
        inds <- which(is.finite(ub))
        m1 <- length(inds)
        B <- matrix(0, nrow = m1, ncol = n)
        for (k in 1:m1) B[k, inds[k]] <- 1
        A <- rbind(A, B)
        b <- c(b, ub[inds])
    }

    # lower bounds
    if (!is.null(lb)) {
        inds <- which(lb > 0)   # all(lb >= 0)
        m2 <- length(inds)
        B <- matrix(0, nrow = m2, ncol = n)
        for (k in 1:m2) B[k, inds[k]] <- -1
        A <- rbind(A, B)
        b <- c(b, -lb[inds])
    }

    # slack variables and base vector
    m <- nrow(A)            # m <- m_in + m1 + m2
    A <- cbind(A, diag(m))
    cc <- c(cc, rep(0, m))
    v <- c(rep(0, n), b)
    I <- (n+1):(n+m)

    return(list(A = A, b = b, cc = cc, x0 = v, I0 = I))
}


.lp.solve <- function(cc, A, b, x, I, maxiter = maxiter) {
    # Linear program  A x = b, minimize cc * x.
    # x is assumed to be a base vector with index set I.
    stopifnot(is.numeric(A), is.numeric(b), is.numeric(x))

    if (!is.matrix(A))
        stop("Argument 'A' must be a numeric matrix.")
    m <- nrow(A); n <- ncol(A)
    if (length(b) != m || length(x) != n)
        stop("One of arguments 'x' or 'b' does not have correct length.")

    # cat("Start value:", cc %*% x, "\n")
    # check I: must be a strict subset of {1, ..., n}

    iter <- 0
    errno <- 0
    while (iter < maxiter) {
        iter <- iter + 1
        J <- setdiff(1:n, I)
        B <- A[, I]; N <- A[, J]
        x_I <- x[I]; x_J <- x[J]
        c_I <- cc[I]; c_J <- cc[J]
        
        # Transform base vector
        y <- solve(t(B), c_I)
        u <- numeric(n)
        for (j in J)
	        u[j] <- cc[j] - t(A[, j]) %*% y
        
        # Case 1: x is a solution of the LP
        if (all(u >= 0)) {
            # cat("Solution found!\n")
            errno <- 1
            break
        }
        
        # Case 2: r in J with u[j] < 0
        r <- which(u < 0)
        if (length(r) > 1) r <- r[sample(1:length(r), 1)]
        
        dd <- solve(B, A[, r])
        d <- numeric(n)
        d[I] <- dd
        
        # Case 2a: The LP does not have a solution
        if (all(d <= 0)) {
            # cat("LP is infeasible!\n")
            errno <- -2
            break
        }
        
        # Case 2b: i in I with d[i] > 0
        i0 <- which(d > 0)
        t <- min(x[i0]/d[i0])
        s <- which(x[i0]/d[i0] == t); s <- i0[s]
        if (length(s) > 1)
            s <- s[sample(1:length(s), 1)]
        
        z <- numeric(n)
        z[r] <- t
        z[I] <- x_I - t*d[I]

        # Define a new base vector
        x <- z
        is <- which(I == s)
        I[is] <- r
        # cat("New base vector:", cc %*% x, "\n")
    }
    if (iter >= maxiter) {
        errno <- -1
    }
    if (errno < 0) {
        x <- y <- NA
    }

    return(list(x = x, dual.x = y, errno = errno))
}