Linear optimization problem with user-defined cost function

Question

I have a linear optimization problem for which I am looking for a suitable optimization solution that can fulfill my requirements. Here is an explanation of the optimization problem:

There are a number of n variables x1, x2, .., xn relevant for the optimization. They contain a value between 0 and 100, whereby the sum of the variables is always 100%. The variables are used as multiplication coefficient for their respective vector with m entries.

There are different functions for determining the cost value for the variables depending on the optimization goal. One of them is shown below as a pseudo code:

cost_function(x1, x2, .., xn)
       i = 0
       min_value = 0
       **while** i < m
           e = x0 * E0_i + x1 * E1_i + .. + xn * EN_i
           v = (s0 * E0_i + s1 * E1_i + .. + sn * EN_i) / e
           minimize_value += min(-C_i, e) * (P_i - v)
       **return** minimize_value

where the bold characters represent vectors and the other scalar values.

The diagram shows the state space with three variables x1, x2 and x3, where darker is better:

Currently, I use Simulated Annealing for this optimization problem, but the state s' is only accepted as new state if it is better than s, because it would be counterproductive to have a probability to take over worse states in a linear optimization problem. From a scientific point of view, the use of Simulated Annealing is probably not appropriate for the application.

Would linear programming still be recommended in this case? Or is there an alternative to Simulated Annealing for linear optimization problems?

---

EDIT:

A minimum working example of the problem can be found here

The result of the variables should be: [0.17023975 0.19532928 0.64296507]

Answer 1

First, the problem is not a linear optimization problem, at least not for the objective function shown (which is nonlinear due the conditional portion in lines 10-13 and particularly the division by E_ges_i in line 13. Simulated annealing might be fine as a heuristic approach, but given the nonlinear objective, only accepting improving steps might or might not be appropriate. If the objective is unimodal, I think skipping steps that make the objective worse might be okay. Then again, if the heat map is unimodal, constrained gradient descent might work just as well or better.

Sticking to the objective function shown, another possibility would be to treat the problem as a mixed integer linear program, using a piecewise linear approximation to the objective function. Since an approximation is involved, you might want to try gradient descent or other local search from the model's solution to see if you can improve it, or alternatively refine the approximation near the model solution and solve the revised model (ad nauseum).

ADDENDUM: I just hacked a little R code, using the objective function posted on GitHub. The algorithm I used is the Nelder-Mead polytope algorithm, as described in 1. The Nelder-Mead algorithm, published in 1965, was a modification of an earlier algorithm (by Spendley, Hext and Himsworth in 1962, according to 1) and has likely been improved upon in turn since the '60s. The attraction of the Nelder-Mead algorithm (or its relatives) here is that it starts with a simplex as the search region and progressively shrinks to smaller and smaller simplices, all of which are subsets of the original. So if we use N-M starting with the unit simplex, the constraints are always satisfied.

I did not put in much care in the coding, did not code all the possible adjustments (such as random restarts), and picked the expansion, contraction and reflection parameters out of thin air. That said, it converged (or at least stopped) after seven iterations on the test problem, with a solution of (0.0126960, 0.2461336, 0.7411704) and an objective value of -4675913 -- not optimal, but I think not too shabby (particularly since I did not code restarts or other more recent tweaks to Nelder-Mead). It's not hard to code, and it does not require derivatives, just function evaluations.

1 P. Gill, W. Murray and M. Wright, Practical Optimization (Academic Press, 1981).

ADDENDUM 2 I updated my R code to use the simplified version of the cost function from the GitHub repository (which is a bit better behaved when an argument is zero). I also switched from the Nelder-Mead algorithm as present in Gill, Murray and Wright to the version on the Wikipedia page, and adjusted the parameter values to the ones they recommend. I have to qualify my earlier comment about Nelder-Mead automatically maintaining feasibility. The condition that the weights sum to 1 is automatically maintained. Nonnegativity of the weights requires occasionally shrink a proposed step, which is easily handled.

With the modified code and the parameter values from the Wikipedia page, I get a final solution of (0, 0.2885719, 0.7114281) with an objective value of -4,683,095.

ADDENDUM 3 Hopefully this will be my final addendum. :-) I also tried a genetic algorithm (in R). The solution from a GA is inherently random, but with the parameters and random seed that I chose I got a final solution of (0.001182725, 0.2869361, 0.7118812) with objective value -4,683,175, which is slightly better than both what I got with Nelder-Mead and what LocalSolver reported.

My R code for both Nelder-Mead and the GA (using the genalg R library) are available in an R notebook.

Answer 2

If you want to implement an algorithm on your own, then we suggest a randomized, derivative-free search, even simpler than a Nelder-Mead approach. Given a feasible solution (respecting the sum equal to 1), move randomly the values of the variables by an epsilon while maintaining the constraint feasible. If the solution is better, keep it; otherwise, throw it away. Start with this simple approach. To go further, tune how you choose the epsilons to move, accept less good solutions along the search to diversify (as done in simulated annealing), and restart the search.

Hexaly (ex LocalSolver), our global optimization solver, combines several optimization techniques under the hood. Here, the above essentially allows Hexaly to perform very well on your problem. Thanks to the small number of dimensions (only 3 variables), there is no need to use derivatives (even approximated) to guide and speed up the search. In the same way, there is no need for surrogate modeling of the cost function because this one is extremely fast to evaluate (about 10,000 calls per second).

Disclaimer: Hexaly is commercial software. You can try it for free for 1 month. In addition, Hexaly is free for basic research and teaching.

Please find below the results obtained by Hexaly using your cost function as external function:

function model() {
    X[0..2] <- float(0,1);
    constraint sum[i in 0..2](X[i]) == 1;
    func <- doubleExternalFunction(cost);
    obj <- call(func, X);
    minimize obj;
}

Having declared the cost function, Hexaly solves the problem as is. Here "solve" means that Hexaly will try to find the best feasible solution to the problem. You can also specify lower and upper bounds for the cost function so that Hexaly computes an optimality gap and then possibly proves the optimality of the solution found.

You can write your model using the Hexaly modeling language (namely LSP) or using Python, Java, C#, or C++ APIs. Here is the link to download the LSP file: https://www.hexaly.com/misc/emma.lsp. Having installed Hexaly, you can run it using the command "hexaly emma.lsp" in console. The best solution found by Hexaly after a few seconds on a basic laptop is:

cost = -4683181.09020784, X0 = 0.00106356929433748, X1 = 0.287235884100486, X2 = 0.711702039130129

The sum over the X equals 1.00000149252495, slightly above 1, because Hexaly uses a tolerance to satisfy constraints. If you wish the sum over the X to be surely lower than 1, then you can set "< 1" in the above model instead of "== 1". In this case, you will find the following solution:

cost = -4683175.50600612, X0 = 0.00111513425966878, X1 = 0.286966585180356, X2 = 0.711915927974678

Now, the sum over the X is equal to 0.999997647414703.