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

______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to