Dear All, I am trying to solve the following maximization problem with R: find x(t) (continuous) that maximizes the integral of x(t) with t from 0 to 1, subject to the constraints dx/dt = u, |u| <= 1, x(0) = x(1) = 0. The analytical solution can be obtained easily, but I am trying to understand whether R is able to solve numerically problems like this one. I have tried to find an approximate solution through discretization of the objective function but with no success so far. Thanks in advance, Paul
On 06/01/2008 7:05 PM, Paul Smith wrote:> Dear All, > > I am trying to solve the following maximization problem with R: > > find x(t) (continuous) that maximizes the > > integral of x(t) with t from 0 to 1, > > subject to the constraints > > dx/dt = u, > > |u| <= 1, > > x(0) = x(1) = 0. > > The analytical solution can be obtained easily, but I am trying to > understand whether R is able to solve numerically problems like this > one. I have tried to find an approximate solution through > discretization of the objective function but with no success so far.R doesn't provide any way to do this directly. If you really wanted to do it in R, you'd need to choose some finite dimensional parametrization of u (e.g. as a polynomial or spline, but the constraint on it would make the choice tricky: maybe a linear spline?), then either evaluate the integral analytically or numerically to give your objective function. Then there are some optimizers available, but in my experience they aren't very good on high dimensional problems: so your solution would likely be quite crude. I'd guess you'd be better off in Matlab, Octave, Maple or Mathematica with a problem like this. Duncan Murdoch
This can be discretized to a linear programming problem so you can solve it with the lpSolve package. Suppose we have x0, x1, x2, ..., xn. Our objective (up to a multiple which does not matter) is: Maximize: x1 + ... + xn which is subject to the constraints: -1/n <= x1 - x0 <= 1/n -1/n <= x2 - x1 <= 1/n ... -1/n <= xn - x[n-1] <= 1/n and x0 = xn = 0 On Jan 6, 2008 7:05 PM, Paul Smith <phhs80 at gmail.com> wrote:> Dear All, > > I am trying to solve the following maximization problem with R: > > find x(t) (continuous) that maximizes the > > integral of x(t) with t from 0 to 1, > > subject to the constraints > > dx/dt = u, > > |u| <= 1, > > x(0) = x(1) = 0. > > The analytical solution can be obtained easily, but I am trying to > understand whether R is able to solve numerically problems like this > one. I have tried to find an approximate solution through > discretization of the objective function but with no success so far. > > Thanks in advance, > > Paul > > ______________________________________________ > R-help at 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. >
On Jan 7, 2008 1:32 AM, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:> This can be discretized to a linear programming problem > so you can solve it with the lpSolve package. Suppose > we have x0, x1, x2, ..., xn. Our objective (up to a > multiple which does not matter) is: > > Maximize: x1 + ... + xn > > which is subject to the constraints: > > -1/n <= x1 - x0 <= 1/n > -1/n <= x2 - x1 <= 1/n > ... > -1/n <= xn - x[n-1] <= 1/n > and > x0 = xn = 0 > > > On Jan 6, 2008 7:05 PM, Paul Smith <phhs80 at gmail.com> wrote: > > Dear All, > > > > I am trying to solve the following maximization problem with R: > > > > find x(t) (continuous) that maximizes the > > > > integral of x(t) with t from 0 to 1, > > > > subject to the constraints > > > > dx/dt = u, > > > > |u| <= 1, > > > > x(0) = x(1) = 0. > > > > The analytical solution can be obtained easily, but I am trying to > > understand whether R is able to solve numerically problems like this > > one. I have tried to find an approximate solution through > > discretization of the objective function but with no success so far.Thats is clever, Gabor! But suppose that the objective function is integral of sin( x( t ) ) with t from 0 to 1 and consider the same constraints. Can your method be adapted to get the solution? Paul
On Jan 6, 2008 8:43 PM, Paul Smith <phhs80 at gmail.com> wrote:> > On Jan 7, 2008 1:32 AM, Gabor Grothendieck <ggrothendieck at gmail.com> wrote: > > This can be discretized to a linear programming problem > > so you can solve it with the lpSolve package. Suppose > > we have x0, x1, x2, ..., xn. Our objective (up to a > > multiple which does not matter) is: > > > > Maximize: x1 + ... + xn > > > > which is subject to the constraints: > > > > -1/n <= x1 - x0 <= 1/n > > -1/n <= x2 - x1 <= 1/n > > ... > > -1/n <= xn - x[n-1] <= 1/n > > and > > x0 = xn = 0 > > > > > > On Jan 6, 2008 7:05 PM, Paul Smith <phhs80 at gmail.com> wrote: > > > Dear All, > > > > > > I am trying to solve the following maximization problem with R: > > > > > > find x(t) (continuous) that maximizes the > > > > > > integral of x(t) with t from 0 to 1, > > > > > > subject to the constraints > > > > > > dx/dt = u, > > > > > > |u| <= 1, > > > > > > x(0) = x(1) = 0. > > > > > > The analytical solution can be obtained easily, but I am trying to > > > understand whether R is able to solve numerically problems like this > > > one. I have tried to find an approximate solution through > > > discretization of the objective function but with no success so far. > > Thats is clever, Gabor! But suppose that the objective function is > > integral of sin( x( t ) ) with t from 0 to 1 > > and consider the same constraints. Can your method be adapted to get > the solution?If a linear approx is sufficient then yes; otherwise, no. For example, if x can be constrained to be small then its roughly true that sin(x) = x and you are back to the original problem.
Hi Paul, Your problem statement does not make much sense to me. You say that an analytical solution can be found easily. I don't see how. This is a variational calculus type problem, where you maximize a functional. Your constraint dx/dt=u(t) means that there exists a solution (the anti-derivative of u) that is unique up to an arbitrary constant. However, a solution may not even exist since you are imposing two conditions on it: x(0) = x(1) = 0. If your solution satisfies both conditions, then it certainly is unique, and it is the x(t) that maximizes integral. Ravi. ---------------------------------------------------------------------------- ------- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: rvaradhan at jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html ---------------------------------------------------------------------------- -------- -----Original Message----- From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Paul Smith Sent: Sunday, January 06, 2008 7:06 PM To: r-help Subject: [R] Can R solve this optimization problem? Dear All, I am trying to solve the following maximization problem with R: find x(t) (continuous) that maximizes the integral of x(t) with t from 0 to 1, subject to the constraints dx/dt = u, |u| <= 1, x(0) = x(1) = 0. The analytical solution can be obtained easily, but I am trying to understand whether R is able to solve numerically problems like this one. I have tried to find an approximate solution through discretization of the objective function but with no success so far. Thanks in advance, Paul ______________________________________________ R-help at 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.
On Jan 8, 2008 2:58 AM, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:> > Thanks again, Gabor. Apparently, something is missing in your > > approach, as I tried to apply it to the original problem (without sin) > > and the result seems not being correct:: > > > > f <- function(z) { > > a <- sum(-(999:1)*z) > > return(a) > > } > > > > result <- optim(rep(0,999), f, > > NULL,lower=rep(-1/1000,999),upper=rep(1/1000,999),method="L-BFGS-B") > > b <- result$par > > > > x <- 0 > > > > for (i in 1:999) { > > x[i+1] <- b[i] + x[i] > > } > > > > The vector x contains the solution, but it does not correspond to the > > analytical solution. The analytical solution is > > > > x(t) = t, if t <= 1/2 and > > > > x(t) = 1 - t, if t >= 1/2. > > > > Am I missing something? > > > > No, but I was -- the x[n] = 0 constraint. Add it through > a penalty term and try this: > > theta <- 1000 > n <- 100 > x <- rep(1/n, n) > f <- function(x) sum(n:1 * x) - theta * sum(x)^2 > optim(x, f, lower = rep(-1, n), upper = rep(1, n), method = "L-BFGS-B", > control = list(maxit = 1000, fnscale = -1, lmm = 25)) > plot(cumsum(res$par))Thanks, Gabor. Inspired by your code, I have written the following: theta <- 50000000 n <- 1000 x <- rep(1/n, n) f <- function(x) sum(n:1 * x) - theta * sum(x)^2 res <- optim(x, f, lower = rep(-1/n, n), upper = rep(1/n, n), method "L-BFGS-B",control = list(maxit = 2000, fnscale = -1)) b <- res$par x <- 0 for (i in 1:n) { x[i+1] <- b[i] + x[i] } plot(seq(0,1,by=1/n),x,type="l") to solve the same problem. It works fine unless one replaces theta <- 50000000 by theta <- 1e100 According to a book that I consulted meanwhile, the larger the value of theta the better the solution should be. Why does not it happen with my example? Any clues? Paul
If its too large the the primary portion of the objective basically just looks like roundoff error. By the way there is also rdonlp2 which can handle constraints. Not on CRAN -- try google. On Jan 9, 2008 6:45 AM, Paul Smith <phhs80 at gmail.com> wrote:> > On Jan 8, 2008 2:58 AM, Gabor Grothendieck <ggrothendieck at gmail.com> wrote: > > > Thanks again, Gabor. Apparently, something is missing in your > > > approach, as I tried to apply it to the original problem (without sin) > > > and the result seems not being correct:: > > > > > > f <- function(z) { > > > a <- sum(-(999:1)*z) > > > return(a) > > > } > > > > > > result <- optim(rep(0,999), f, > > > NULL,lower=rep(-1/1000,999),upper=rep(1/1000,999),method="L-BFGS-B") > > > b <- result$par > > > > > > x <- 0 > > > > > > for (i in 1:999) { > > > x[i+1] <- b[i] + x[i] > > > } > > > > > > The vector x contains the solution, but it does not correspond to the > > > analytical solution. The analytical solution is > > > > > > x(t) = t, if t <= 1/2 and > > > > > > x(t) = 1 - t, if t >= 1/2. > > > > > > Am I missing something? > > > > > > > No, but I was -- the x[n] = 0 constraint. Add it through > > a penalty term and try this: > > > > theta <- 1000 > > n <- 100 > > x <- rep(1/n, n) > > f <- function(x) sum(n:1 * x) - theta * sum(x)^2 > > optim(x, f, lower = rep(-1, n), upper = rep(1, n), method = "L-BFGS-B", > > control = list(maxit = 1000, fnscale = -1, lmm = 25)) > > plot(cumsum(res$par)) > > Thanks, Gabor. Inspired by your code, I have written the following: > > theta <- 50000000 > n <- 1000 > x <- rep(1/n, n) > f <- function(x) sum(n:1 * x) - theta * sum(x)^2 > res <- optim(x, f, lower = rep(-1/n, n), upper = rep(1/n, n), method > "L-BFGS-B",control = list(maxit = 2000, fnscale = -1)) > > b <- res$par > > x <- 0 > > for (i in 1:n) { > x[i+1] <- b[i] + x[i] > } > > plot(seq(0,1,by=1/n),x,type="l") > > to solve the same problem. It works fine unless one replaces > > theta <- 50000000 > > by > > theta <- 1e100 > > According to a book that I consulted meanwhile, the larger the value > of theta the better the solution should be. Why does not it happen > with my example? Any clues? > > Paul > > > ______________________________________________ > R-help at 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. >