# Linear objective function with non-linear constraints

I would like to choose a set of $$\beta_j$$s that maximizes a simple linear objective function of the type

$$\underset{\beta_j}{\operatorname{max}}\sum_{j=1}^{J}X_j\beta_j \\$$

subject to the following constraints $$\sum_{j=1}^{J}C_j(\beta_j)\beta_j \le M \\ \beta_j \in \Omega \\$$

here $$C_j(\beta_j)$$ can be thought as a marginal cost function that changes with the chosen $$\beta_j$$. $$\beta_j$$ can only be from a set of pre-selected set of integers $$\Omega$$. $$M$$ is some budget constraint.

I don't know the functional form of $$C_j(\beta_j)$$ but can simulate $$C_j$$ for each $$j$$ and each possible $$\beta_j$$.

I am having trouble understanding how to optimize this problem efficiently. Can someone give any direction on how this can be solved in R or Python?

• Given that you don't have that many possible values for $\beta_j$, could you simply replace the $\beta_j$ with binary $x_{ij}$ equalling 1 iff $\beta_j$ equals the $i$'th value in $\Omega$ – Sune Jul 15 at 18:06
• Let $\omega_i$, $i\in I$, be the different values for the $\beta$ variables. Then, if I understand you correctly, you can precompute $C_j(\omega_i)=\gamma_{ij}$ for all &$i$ and $j$. Now, replace the $\beta_j$-variables with $\sum_{i\in I} \omega_ix_{ij}$. You need constraints of the form $\sum_{i\in I} x_{ij} =1$ for all $j$ and your budget constraint should be $\sum_{j=1}^J\sum_{i\in I} \gamma_{ij} \omega_ix_{ij} \leq M$ (if I'm not mistaken) – Sune Jul 15 at 19:50
• I.e., you would end up with a MILP, for which there are many solvers. – Mark L. Stone Jul 16 at 3:52
• There are many options for formulating and solving MILPs in R. One of the easiest might be using CVXR cran.r-project.org/web/packages/CVXR/index.html and cvxr.rbind.io However, this might not provide access to all the low-level control functionality available in "native" interfaces. CVXPY cvxpy.org is similar, but for Python rather than R. – Mark L. Stone Jul 16 at 14:00
• Yes, but I have no specific experience in or knowledge of it. What does the company not like about installing CVXR (note that this is not CVX, which is a package for MATLAB, confusingly based at the domian cvxr.com , standing for CVX Research)) . Actually CVX is the granddaddy of them all, and CVXPY and CVXR were later developed by other students of the same Prof (Stephen Boyd) who was an initial co-developer of CVX along with his student) – Mark L. Stone Jul 16 at 14:50

Given that my comment to the question has been referred in the other existing answer, I will add it as an answer on its own. The premise for the answer is that the $$C_j(y)$$ function values can be precomputed for all values of $$y\in\Omega$$. The basic idea is to utilize that $$\Omega$$ only contains a relatively small number of values in order to convert the problem into a binary linear program solvable by many commercial as well as free solvers.

To that end, let $$\omega_i$$, $$i\in I$$, be the different values in $$\Omega$$. Then, for each $$i \in I$$ and $$j\in J$$ compute the values $$C_j(\omega_i):=\gamma_{ij}$$. Next, introduce binary variables $$z_{ij}$$ equaling 1 iff $$\beta_j$$ takes the value $$\omega_i$$. We can then replace the variables $$\beta_j$$ with the sum $$\sum_{i\in I}\omega_iz_{ij}$$. The original problem can the be stated as \begin{align} \max& \sum_{j\in J}X_j \sum_{i \in I}\omega_iz_{ij}\\ \text{s.t.}:&\sum_{i\in I} z_{ij} = 1,&& \forall j\in J\\ & \sum_{j\in J} \sum_{i\in I} \gamma_{ij}\omega_iz_{ij}\leq M,\\ & z_{ij}\in \{0,1\},&&\forall i\in I,j \in J \end{align} I could imagine, but haven’t tested it, that many solvers can handle this somewhat simple MILP efficiently.

Since you don't know the functional form, you can use Pypopt, the Python wrapper around Ipopt. Ipopt supports callbacks, which means you can provide functions for the solver to evaluate in real-time to get values & derivatives.

Another way would be to use any of the genetic/evolutionary algorithms in Scipy.

If you have values in a tabular format, i.e., you don't have a black-box function that can produce $$C(\beta)$$ for any $$\beta$$, the workaround for non-linear optimisation would be to simply interpolate between the closest values that you do have. Ipopt defaults to finite differences if you don't provide derivatives, so as a first-order approach you would only need to do this for the evaluation of the function (not the derivatives).

It's important to know that it's incorrect to solve this directly as an MILP, as your $$C(\beta)$$ would be fixed, rather than updated dynamically as it's supposed to.

If you do want to use an MILP formulation to select values from a table you can, but with a few caveats:

• you would lose the derivative-based information
• The number of new binaries won't scale well for dense tables
• The formulation can be quite challenging
• There's a decent chance you will need a commercial linear solver

Thus, the best all-round (and free) option in my opinion would be callbacks through Ipopt.

• I thought @Sune was suggesting pre-computing for all possible values. Therefore, the "dynamic updating" would be accomplished via (optimizer choosing the value of) the binary variables. I.e.,, the resulting MILP would have what I think you are referring to as "dynamic updating". – Mark L. Stone Jul 16 at 13:50
• @MarkL.Stone Yeah that could work as long as they can somehow feed the current values of the function into the solver as callbacks - not sure if linear solvers support this though. If done through a table formulation (which is doable) it might be a bit too much for linear solvers to handle. – Nikos Kazazakis Jul 16 at 15:11
• @NikosKazazakis Thanks for the solution. The main issue that I face with this is that I don't have access to the simulation that provides the $C_j$s. I have been a matrix of $C_j$s for all $j$ and $\Omega$. – FightMilk Jul 16 at 15:45
• @FightMilk In practice you could approximate that through the callbacks, as in if you don't have an exact $C$ for a certain $\beta$ you can interpolate between the two closest values. As Sune suggested you could use an MILP formulation if that's viable size/complexity-wise. In that case however, you would lose the derivative-based information that the nonlinear solver would otherwise have. – Nikos Kazazakis Jul 16 at 15:56
• Thank you and that makes a lot of sense. I have very new to MILP and optimization in general. I will do the necessary reading and on Pyopt package and get back! – FightMilk Jul 16 at 16:08