Liu Evans, Gareth
2011-Oct-03 15:29 UTC
[R] minimisation problem, two setups (nonlinear with equality constraints/linear programming with mixed constraints)
Dear All,
Thank you for the replies to my first thread here:
http://r.789695.n4.nabble.com/global-optimisation-with-inequality-constraints-td3799258.html.
So far the best result is achieved via a penalised objective function. This was
suggested by someone on this list privately. I am still looking into some of
the options mentioned in the original thread, but I have been advised that there
may be better ways if I present the actual problem with a reproducible example.
In principle the problem can be solved by linear programming, so I include code
for my attempt at this via RGLPK. It says that there is no feasible solution,
but the solution is known analytically in the case below.
Here is the precise problem:
Minimise, over 100?1 real vectors v,
Max_i(|v_i|) such that X'v=e_2,
where X is a given 100?2 matrix and e_2 =(0,1)'. The v_i are the elements
of v.
I have put the actual X matrix at the end of this post, along with a feasible
starting value for v.
The correct minimum is 0.01287957, obtained with v_i=0.01287957 for i<=50 and
v_i = 0.01287957 for i>=51.
Here is the code for the penalised objective function approach, adapted very
slightly from what someone on this list sent me:
......................................................................
X <- # See end of this message for the X data
x1 <- X[, 1]
x2 <- X[, 2]
fun <- function(q) {
mu <- 0.1
max(abs(v)) + (sum(v*x1)^2 + (1-sum(x2*v))^2)/(2*mu)
}
vstart <- # feasible starting value. See end of this post.
sol <- optim(vstart, fun, method="L-BFGS-B", lower=rep(-1,
100), upper=rep(1,100))
max(abs(sol$par))
.........................................................................
This gets quite near, around 0.015-0.016 for me by varying mu.
Alternatively the problem can be set up as a linear programming task as follows:
Minimise, over 100?1 real vectors v, and over scalars q >= 0,
q
such that, for i=1,...,100
v_i<=q
v_i>=-q
X'v=e_2
Here is my RGLPK code:
.........................................................................
X <- # See end of this message for the X data
XROWS <- 100
XCOLS <- 2
e_2=rep(0,times=XCOLS)
e2[2]<- 1
obj <- c(rep(0,XROWS),1) # coefficients on v_1, . . . , v_100, q.
mat <- matrix(rbind(cbind(diag(XROWS), rep(-1,XROWS)), cbind(diag(XROWS),
rep(1,XROWS)), cbind(t(X), rep(0,XCOLS)), cbind(t(rep(0,XROWS)), 1)),
nrow=2*XROWS+XCOLS+1)
dir <- c(rep("<=", XROWS), rep(">=", XROWS),
rep("==", XCOLS), ">=")
rhs <- c(rep(0, 2*XROWS), e_2, 0)
sol <- Rglpk_solve_LP(obj, mat, dir, rhs, types = NULL, max = FALSE,
bounds = c(-5,5), verbose = TRUE)
...........................................................................
The output is
"
GLPK Simplex Optimizer, v4.42
203 rows, 101 columns, 601 non-zeros
0: obj = 0.000000000e+000 infeas = 1.000e+000 (2)
4: obj = 0.000000000e+000 infeas = 1.000e+000 (1)
PROBLEM HAS NO FEASIBLE SOLUTION
"
I have also tried setting the problem up with a small interval around the
equality constraints rather than having strict equalities, but could not get the
correct solution this way either.
Maybe I am making an error with RGLPK - I have been told it should work for a
problem of this size and much larger. I have also tried DEoptim, IpSolve and
ConstrOptim.
Regards,
Gareth
...............................................................................
Values of X and vstart below
...............................................................................
vstart
-0.025251183
-0.022301089
-0.020429759
-0.01902228
-0.017877586
-0.016903415
-0.016049376
-0.015284788
-0.014589517
-0.0139496
-0.01335494
-0.012797983
-0.012272922
-0.011775194
-0.011301139
-0.010847778
-0.010412649
-0.009993692
-0.009589163
-0.009197573
-0.00881764
-0.008448248
-0.008088421
-0.007737298
-0.007394116
-0.007058194
-0.006728922
-0.006405748
-0.006088173
-0.005775744
-0.005468043
-0.005164689
-0.00486533
-0.004569638
-0.00427731
-0.003988063
-0.003701629
-0.003417759
-0.003136215
-0.002856772
-0.002579217
-0.002303343
-0.002028954
-0.00175586
-0.001483876
-0.001212825
-0.00094253
-0.00067282
-0.000403526
-0.000134481
0.000134481
0.000403526
0.00067282
0.00094253
0.001212825
0.001483876
0.00175586
0.002028954
0.002303343
0.002579217
0.002856772
0.003136215
0.003417759
0.003701629
0.003988063
0.00427731
0.004569638
0.00486533
0.005164689
0.005468043
0.005775744
0.006088173
0.006405748
0.006728922
0.007058194
0.007394116
0.007737298
0.008088421
0.008448248
0.00881764
0.009197573
0.009589163
0.009993692
0.010412649
0.010847778
0.011301139
0.011775194
0.012272922
0.012797983
0.01335494
0.0139496
0.014589517
0.015284788
0.016049376
0.016903415
0.017877586
0.01902228
0.020429759
0.022301089
0.025251183
.....................................................................................................................................
X
1 -2.330078923
1 -2.057855981
1 -1.885177032
1 -1.755300501
1 -1.649672679
1 -1.559779992
1 -1.480972651
1 -1.410419531
1 -1.346262665
1 -1.287213733
1 -1.232340861
1 -1.180947041
1 -1.13249653
1 -1.086568115
1 -1.042824239
1 -1.000989917
1 -0.960837931
1 -0.922178178
1 -0.884849841
1 -0.848715527
1 -0.813656808
1 -0.779570774
1 -0.746367337
1 -0.713967098
1 -0.682299633
1 -0.651302112
1 -0.62091817
1 -0.591096977
1 -0.561792466
1 -0.532962693
1 -0.504569287
1 -0.476576998
1 -0.448953298
1 -0.421668052
1 -0.39469322
1 -0.368002611
1 -0.341571661
1 -0.315377237
1 -0.289397474
1 -0.263611615
1 -0.237999879
1 -0.212543342
1 -0.187223821
1 -0.162023779
1 -0.136926226
1 -0.111914639
1 -0.08697288
1 -0.062085116
1 -0.037235755
1 -0.012409369
1 0.012409369
1 0.037235755
1 0.062085116
1 0.08697288
1 0.111914639
1 0.136926226
1 0.162023779
1 0.187223821
1 0.212543342
1 0.237999879
1 0.263611615
1 0.289397474
1 0.315377237
1 0.341571661
1 0.368002611
1 0.39469322
1 0.421668052
1 0.448953298
1 0.476576998
1 0.504569287
1 0.532962693
1 0.561792466
1 0.591096977
1 0.62091817
1 0.651302112
1 0.682299633
1 0.713967098
1 0.746367337
1 0.779570774
1 0.813656808
1 0.848715527
1 0.884849841
1 0.922178178
1 0.960837931
1 1.000989917
1 1.042824239
1 1.086568115
1 1.13249653
1 1.180947041
1 1.232340861
1 1.287213733
1 1.346262665
1 1.410419531
1 1.480972651
1 1.559779992
1 1.649672679
1 1.755300501
1 1.885177032
1 2.057855981
1 2.330078923
Liu Evans, Gareth
2011-Oct-07 13:25 UTC
[R] Linear programming problem, RGPLK - "no feasible solution"
Dear All,
RGLPK says that there is no feasible solution but I think there should be. In a
more general setup with more variables, where we know the feasible solution
analytically, it also says there is no feasible solution. I have tried
different values. except the 0s and 1s must stay as they are. The result has
always been the same so far, i.e. no feasible solution.
Here is my code:
obj <- c(rep(0,3),1)
col1 <-c(1,0,0,1,0,0,1,-2.330078923,0)
col2 <-c(0,1,0,0,1,0,1,-2.057855981,0)
col3 <-c(0,0,1,0,0,1,1,-1.885177032,0)
col4 <-c(-1,-1,-1,1,1,1,0,0,1)
mat <- cbind(col1, col2, col3, col4)
dir <- c(rep("<=", 3), rep(">=", 3),
rep("==", 2), ">=")
rhs <- c(rep(0, 6), ej, 0)
sol <- Rglpk_solve_LP(obj, mat, dir, rhs, types = NULL, max = FALSE,
bounds = c(-100,100), verbose = TRUE)
Regards,
Gareth