2015-05-27 14 views
6

Sto cercando di trovare il minimo locale di una funzione ei parametri hanno una somma fissa. Ad esempio,Ottimizzazione R con vincoli di uguaglianza e disuguaglianza

Fx = 10 - 5x1 + 2x2 - x3

e le condizioni sono le seguenti,

x1 + x2 + x3 = 15

(x1, x2, x3)> = 0

Dove la somma di x1, x2 e x3 ha un valore noto e sono tutti maggiori di zero. In R, sarebbe simile a questa,

Fx = function(x) {10 - (5*x[1] + 2*x[2] + x[3])} 
opt = optim(c(1,1,1), Fx, method = "L-BFGS-B", lower=c(0,0,0), upper=c(15,15,15)) 

Ho anche cercato di usare le disuguaglianze con constrOptim per forzare la somma da fissare. Penso ancora che questo possa essere un lavoro plausibile, ma non sono riuscito a farlo funzionare. Questo è un esempio semplificato del problema reale, ma qualsiasi aiuto sarebbe molto apprezzato.

risposta

6

In questa occasione optim non funzionerà ovviamente perché si dispone di vincoli di uguaglianza. constrOptim non funzionerà per lo stesso motivo (ho provato a convertire l'uguaglianza in due disuguaglianze, cioè maggiore e minore di 15, ma questo non ha funzionato con constrOptim).

Tuttavia, esiste un pacchetto dedicato a questo tipo di problema e cioè Rsolnp.

si utilizza la seguente modo:

#specify your function 
opt_func <- function(x) { 
    10 - 5*x[1] + 2 * x[2] - x[3] 
} 

#specify the equality function. The number 15 (to which the function is equal) 
#is specified as an additional argument 
equal <- function(x) { 
    x[1] + x[2] + x[3] 
} 

#the optimiser - minimises by default 
solnp(c(5,5,5), #starting values (random - obviously need to be positive and sum to 15) 
     opt_func, #function to optimise 
     eqfun=equal, #equality function 
     eqB=15, #the equality constraint 
     LB=c(0,0,0), #lower bound for parameters i.e. greater than zero 
     UB=c(100,100,100)) #upper bound for parameters (I just chose 100 randomly) 

uscita:

> solnp(c(5,5,5), 
+  opt_func, 
+  eqfun=equal, 
+  eqB=15, 
+  LB=c(0,0,0), 
+  UB=c(100,100,100)) 

Iter: 1 fn: -65.0000  Pars: 14.99999993134 0.00000002235 0.00000004632 
Iter: 2 fn: -65.0000  Pars: 14.999999973563 0.000000005745 0.000000020692 
solnp--> Completed in 2 iterations 
$pars 
[1] 1.500000e+01 5.745236e-09 2.069192e-08 

$convergence 
[1] 0 

$values 
[1] -10 -65 -65 

$lagrange 
    [,1] 
[1,] -5 

$hessian 
      [,1]  [,2]  [,3] 
[1,] 121313076 121313076 121313076 
[2,] 121313076 121313076 121313076 
[3,] 121313076 121313076 121313076 

$ineqx0 
NULL 

$nfuneval 
[1] 126 

$outer.iter 
[1] 2 

$elapsed 
Time difference of 0.1770101 secs 

$vscale 
[1] 6.5e+01 1.0e-08 1.0e+00 1.0e+00 1.0e+00 

Così i valori ottimali risultanti sono:

$pars 
[1] 1.500000e+01 5.745236e-09 2.069192e-08 

che significa che il primo parametro è 15 e il resto zero e zero. Questo è effettivamente il minimo globale nella tua funzione dato che x2 sta aggiungendo alla funzione e 5 * x1 ha un'influenza molto più grande (negativa) di x3 sul risultato. La scelta di 15, 0, 0 è la soluzione e il minimo globale per la funzione in base ai vincoli.

La funzione ha funzionato benissimo!

4

Questo è in realtà un problema di programmazione lineare, quindi un approccio naturale sarebbe utilizzare un solver di programmazione lineare come il pacchetto lpSolve. È necessario fornire una funzione obiettivo e una matrice dei vincoli e il risolutore farà il resto:

library(lpSolve) 
mod <- lp("min", c(-5, 2, -1), matrix(c(1, 1, 1), nrow=1), "=", 15) 

Quindi è possibile accedere alla soluzione ottimale e il valore obiettivo (aggiungendo il termine costante 10, che non è previsto al risolutore):

mod$solution 
# [1] 15 0 0 
mod$objval + 10 
# [1] -65 

un risolutore di programmazione lineare dovrebbe essere un buon affare più veloce di un risolutore generale ottimizzazione lineare e non dovrebbero avere problemi restituendo la soluzione ottimale esatta (invece di un punto vicina che può essere soggetta ad errori di arrotondamento).

+1

Bello! Quando l'OP dice: "Questo è un esempio semplificato del problema reale" mi fa pensare che il problema reale potrebbe essere non lineare. Quindi, per sicurezza, ho suggerito un metodo non lineare (che funziona comunque anche se è più lento). Fornire il gradiente (semplice per questo caso) rende ancora più veloce se la velocità è un problema. Ad ogni modo, non intendo male, è davvero positivo che tu abbia aggiunto questa risposta, sicuramente utile. – LyzandeR