A simple penalty barrier function is formed which is then evaluated at randomly
sampled points based on the upper and lower parameter bounds
(when `eval.type`

= 2), else the objective function directly for values not
violating any inequality constraints (when `eval.type`

= 1). The sampled
points can be generated from the uniform, normal or truncated normal
distributions.

1 2 3 4 |

`pars` |
The starting parameter vector. This is not required unless the fixed option is also used. |

`fixed` |
The numeric index which indicates those parameters which should stay fixed instead of being randomly generated. |

`fun` |
The main function which takes as first argument the parameter vector and returns a single value. |

`eqfun` |
(Optional) The equality constraint function returning the vector of evaluated equality constraints. |

`eqB` |
(Optional) The equality constraints. |

`ineqfun` |
(Optional) The inequality constraint function returning the vector of evaluated inequality constraints. |

`ineqLB` |
(Optional) The lower bound of the inequality constraints. |

`ineqUB` |
(Optional) The upper bound of the inequality constraints. |

`LB` |
The lower bound on the parameters. This is not optional in this function. |

`UB` |
The upper bound on the parameters. This is not optional in this function. |

`distr` |
A numeric vector of length equal to the number of parameters, indicating the choice of distribution to use for the random parameter generation. Choices are uniform (1), truncated normal (2), and normal (3). |

`distr.opt` |
If any choice in |

`bestN` |
The best N (less than or equal to n.sim) set of parameters to return. |

`n.sim` |
The number of random parameter sets to generate. |

`cluster` |
If you want to make use of parallel functionality, initialize and pass a cluster object from the parallel package (see details), and remember to terminate it! |

`rseed` |
(Optional) A seed to initiate the random number generator, else system time will be used. |

`eval.type` |
Either 1 (default) for the direction evaluation of the function (excluding inequality constraint violations) or 2 for the penalty barrier method. |

`trace` |
(logical) Whether to display the progress of the function evaluation. |

`...` |
(Optional) Additional parameters passed to the main, equality or inequality functions |

Given a set of lower and upper bounds, the function generates, for those
parameters not set as fixed, random values from one of the 3 chosen
distributions. For simple functions with only inequality constraints, the direct
method (`eval.type`

= 1) might work better. For more complex setups with
both equality and inequality constraints the penalty barrier method
(`eval.type`

= 2)might be a better choice.

A matrix of dimension bestN x (no.parameters + 1). The last column is the evaluated function value.

The choice of which distribution to use for randomly sampling the parameter
space should be driven by the user's knowledge of the problem and confidence or
lack thereof of the parameter distribution. The uniform distribution indicates a
lack of confidence in the location or dispersion of the parameter, while the
truncated normal indicates a more confident choice in both the location and
dispersion. On the other hand, the normal indicates perhaps a lack
of knowledge in the upper or lower bounds, but some confidence in the location
and dispersion of the parameter. In using choices (2) and (3) for `distr`

,
the `distr.opt`

list must be supplied with `mean`

and `sd`

as
subcomponents for those parameters not using the uniform.

Alexios Ghalanos and Stefan Theussl

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 | ```
## Not run:
library(Rsolnp)
library(parallel)
# Windows
cl = makePSOCKcluster(2)
# Linux:
# makeForkCluster(nnodes = getOption("mc.cores", 2L), ...)
gofn = function(dat, n)
{
x = dat[1:n]
y = dat[(n+1):(2*n)]
z = dat[(2*n+1):(3*n)]
ii = matrix(1:n, ncol = n, nrow = n, byrow = TRUE)
jj = matrix(1:n, ncol = n, nrow = n)
ij = which(ii<jj, arr.ind = TRUE)
i = ij[,1]
j = ij[,2]
# Coulomb potential
potential = sum(1.0/sqrt((x[i]-x[j])^2 + (y[i]-y[j])^2 + (z[i]-z[j])^2))
potential
}
goeqfn = function(dat, n)
{
x = dat[1:n]
y = dat[(n+1):(2*n)]
z = dat[(2*n+1):(3*n)]
apply(cbind(x^2, y^2, z^2), 1, "sum")
}
n = 25
LB = rep(-1, 3*n)
UB = rep( 1, 3*n)
eqB = rep( 1, n)
sp = startpars(pars = NULL, fixed = NULL, fun = gofn , eqfun = goeqfn,
eqB = eqB, ineqfun = NULL, ineqLB = NULL, ineqUB = NULL, LB = LB, UB = UB,
distr = rep(1, length(LB)), distr.opt = list(), n.sim = 2000,
cluster = cl, rseed = 100, bestN = 15, eval.type = 2, n = 25)
#stop cluster
stopCluster(cl)
# the last column is the value of the evaluated function (here it is the barrier
# function since eval.type = 2)
print(round(apply(sp, 2, "mean"), 3))
# remember to remove the last column
ans = solnp(pars=sp[1,-76],fun = gofn , eqfun = goeqfn , eqB = eqB, ineqfun = NULL,
ineqLB = NULL, ineqUB = NULL, LB = LB, UB = UB, n = 25)
# should get a value of around 243.8162
## End(Not run)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.