[R] minimisation problem, two setups (nonlinear with equality constraints/linear programming with mixed constraints)
Liu Evans, Gareth
Gareth.Liu-Evans at liverpool.ac.uk
Mon Oct 3 17:29:02 CEST 2011
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
More information about the R-help
mailing list