Approach for optimization problem with polynomial constraints?

Question

I have a problem where the objective function is linear and constraints have polynomials (in one variable). So, my question is what are the main approaches to this issue? I can construct a small example, just to illustrate it.

$ \max \sum_{i} a_i x_i - \sum_{j} b_j y_j $

$\qquad c_1 x_i + c_2 x_i^2 + c_3 x_i^3 +\ldots + c_k x_i^k = \sum_{j} d_j y_j, \quad \forall i\in N $

$\qquad x_i \geq 0, \quad i\in N $

$\qquad y_j \in \{0,1\}, \quad j\in M $

Answer 1

For small scale problems, simply using a global solver appears to work very well, at least for the data I tried. Here is some YALMIP code (MATLAB Toolbox, developed by me) to solve a small instance using YALMIPs global solver bmibnb. It is solved in a second or so if you have a good MILP solver installed. Similiar with scips global solver

N = 10;
M = 20;
degree = 4;
a = randn(N,1);
b = randn(M,1);
c = rand(M,degree);
d = randn(M,1);

x = sdpvar(N,1);
y = binvar(M,1);
objective = a'*x - b'*y;
Model = x>=0
for i = 1:N
    Model = [Model, [-d'*y c(i,:)]*monolist(x(i),4) == 0];
end

optimize(Model,objective,sdpsettings('solver','bmibnb'))
value(x)
value(y)
%optimize(Model,objective,sdpsettings('solver','scip'))

EDIT: To follow up on your comment, here is a model based on a PWA approximation (using sos2 constructs in cplex, as that speeds up things). Of course, the bound 5 should be chosen more carefully by performing bound propagation etc. Solved in a fraction of a second, but the drawback is of course the lack of an exact solution

xi = linspace(0,5,200)';
Model = x>=0
for i = 1:N
    fi = c(i,1)*xi + c(i,2)*xi.^2 + c(i,3)*xi.^3 + c(i,4)*xi.^4;
    lambda = sdpvar(length(fi),1);
    Model = [Model,sos2(lambda)];
    Model = [Model, x(i) == lambda'*xi, d'*y == lambda'*fi,lambda>=0,sum(lambda)==1]; 
end