# smplx.r simplex algorithm # J F Monahan, July 2009 # # first the code # # phase 1 / phase 2 algorithm p1p2 <- function(A,b,c,m,n) { p <- n+(1:m) # indices of starting basic feasible solution print(dim(A)) Anew <- cbind(A*sign(b),diag(rep(1,m))) # append identity print(dim(Anew)) cnew <- c(rep(0,n),rep(1,m)) # new cost vector feas <- simplex(p,Anew,abs(b),cnew,m,n+m) # solve lp problem print("******donewithphase1****************") print(feas) print("phase1resultsabove") if( (feas[[3]] == 0)*(feas[[4]] == 0 ) ) { # found feasible soln p <- feas[[1]] p1p2 <- simplex(p,A,b,c,m,n) } else { # no feasible soln print("nofeasiblesolution") p1p2 <- list(feas[[1]],feas[[2]],3,0) } } # end of function p1p2 # # simplex algorithm simplex <- function(p,A,b,c,m,n) { # dumb simplex algorithm # p is vector of indices for basic, feasible soln # A matrix, b rhs vector, c costs # m constraints, n variables #p# print(dim(A)) #p# print(p) # initial solution indices B <- A[,p] # matrix of starting solution BiAb <- solve(B,cbind(A,b)) # start loop for( i in (1:100)) { # begin loop #p# print(p) # current basic solution indices pn <- (1:n)[-p] #p# print(pn) # nonbasic indices cp <- c[p] cn <- c[pn] x <- rep(0,n) x[p] <- BiAb[,n+1] #p# print(BiAb[,n+1]) #p# print(sum(cp*BiAb[,n+1])) # current obj fn value # find feasible directions zjmcj <- cp %*% BiAb[,pn] - cn #p# print(zjmcj) flist <- (1:(n-m))*(zjmcj > 1.e-10) #p# print(flist) # potential directions #p print((BiAb[,flist]>1.e-10)) #p print(apply((BiAb[,flist]> 1.e-10)+0,2,sum)) if( sum(flist) == 0 ) { print("done") code <- 0 break } ibig <- which.max(zjmcj) # ibig is in, what is out? #p# print(c(ibig,max(zjmcj))) ibig <- pn[ibig] #p# print(ibig) # what new variable is in? if( max(BiAb[,ibig]) < 1.e-10 ) { print("unbound") code <- 2 break } v <- BiAb[,n+1]/BiAb[,ibig] v[ (BiAb[,ibig]<1.e-10) ] <- 1.e30 lmin <- which.min( v ) # lmin is out ### print(c(lmin,v[lmin])) v <- BiAb[,ibig]/BiAb[lmin,ibig] v[lmin] <- (BiAb[lmin,ibig]-1)/BiAb[lmin,ibig] BiAb <- BiAb - outer(v,BiAb[lmin,]) ### print(BiAb) p[lmin] <- ibig code <- 1 # too many iterations if leave with this } # end of loop x <- rep(0,n) x[p] <- BiAb[,n+1] # return indices, x, return code, obj fn value simplex <- list(p,x,code,sum(c*x)) } # done with simplex # # # call phase 1/ phase2 algorithm # # first problem from Wagner, p103, 119-120 m <- 3 n <- 7 A <- matrix(c(1,7,3,1,5,5,1,3,10,1,2,15,1,0,0,0,1,0,0,0,1),m,n) b <- c(15,120,100) c <- c(-4,-5,-9,-11,0,0,0) A b c prob1 <- p1p2(A,b,c,m,n) # return indices, x, return code, obj fn value prob1 # # second problem from VJ Bowman m <- 6 n <- 10 A <- cbind( matrix(c(1,2,5,0,0,0,4,3,1,0,0,0,5,0,0,1,0,-3,2,0,0,0,1,-4),6,4), diag(6) ) # block identity in cols 5-10 b <- c(7,6,5,4,3,-12) c <- c(-1,-8,-5,-6,rep(0,6)) A b c prob2 <- p1p2(A,b,c,m,n) # return indices, x, return code, obj fn value prob2 # # third problem from Wagner -- unbounded m <- 2 n <- 4 A <- matrix(c(1,-1,-1,1,1,0,0,1),m,n) b <- c(1,1) c <- c(-1,-1,0,0) A b c prob3 <- p1p2(A,b,c,m,n) # return indices, x, return code, obj fn value prob3 # # fourth problem from Wagner -- no solution # same A, c, but negate b b <- -b A b c prob4 <- p1p2(A,b,c,m,n) # return indices, x, return code, obj fn value prob4 # # fifth problem is L1 regression m <- 6 # number of obs n <- 16 # 2*nobs + 2*(cols of X) X <- cbind(rep(1,m),(1:m)) # design matrix A <- cbind(diag(rep(1,m)),diag(rep(-1,m)),X,-X) b <- c(4.5,5.5,6.5,8,10,12) # y's c <- c(rep(1,2*m),0,0,0,0) A b c p <- (1:m) # we have feasible soln, b positive prob5 <- simplex(p,A,b,c,m,n) # return indices, x, return code, obj fn value prob5 # done rm(list=ls()) q()