How to represent the objective function of the Weapon Target Assignment problem in CVXPY?

Question

I am trying to use CVXPY to analyse a problem and the objective function for this problem involves calculating a product and a sum as per the problem description below (taken from a draft paper I am writing):

---

---

However, I cannot figure out how to specify such an objective function in CVXPY given that it uses a product to combine all of the probability values for a particular target. In order to replicate the objective function, I need to multiply all of the elements of each row in a matrix together (i.e. multiplication along an axis). I am hoping to achieve something like the function shown below, which doesn't work because of the use of the numpy function.

objective = cp.Maximize(cp.sum(np.prod(1-cp.exp(cp.multiply(x, cp.log(q))), axis=1)))

Is there any way that I can perform multiplication along an axis in CVXPY? Or, is there an alternative way to represent my objective function?

Answer 1

TLDR: This can be formulated and solved in CXVPY with Mosek as solver, as a Mixed-Integer (generalized) Geometric Programming problem, using CVXPY's Disciplined Geometric Programming (DGP) capability. It can also be formulated and solved in CVX using its gp mode.

Details: This satisfies all of CVXPY's Disciplined Geometric Programming (DGP) rules. The objective function is a sum of monomials, which is a posynomial, and is log-log convex. So this can be formulated as a (generalized) Mixed-Integer Geometric Programming problem. Mosek is the best (only?) solver to use, because it has native exponential cone capability, which is invoked by CVXPY to solve the generalized Geometric Programming problem, and can solve such problems subject to integer constraints on some or all of the variables.

See the CVXPY web page for Disciplined Geometric Programming (DGP),

Answer 2

Now $$ \begin{array}{rcl} t_j & \geq & \prod_{i=1}^m (1-p_{ij})^{x_{ij}} \\ & = & exp(\ln(\prod_{i=1}^m (1-p_{ij})^{x_{ij}})) \\ & = & exp(\sum_{i=1}^m x_{ij} \ln(1-p_{ij})) \\ \end{array} $$ This is the same as saying $$ (t_j;1;\sum_{i=1}^m x_{ij} \ln(1-p_{ij})) $$ has to belong to the primal exponential cone.

So you can easily rewrite your problem as a conic optimization problem involving the exponential cone.

Mosek for sure can solve those problems even if some variables are constrained to be integers. I guess Cvxpy+Mosek should be handle that formulation. Not sure how hard the problem will be given the integer constraints.

Answer 3

There is no need for DGP rules - the problem can be formulated directly as a mixed integer DCP, and you were pretty close to it. I know you said you figured it out, but I'm adding this for future visits to this question (since neither give a CVXPY implementation).

Somewhat in line with your code, I'm setting this up initially as:

$\min \sum_{j=1}^n \prod_{i=1}^m q_{ij}^{x_{ij}}$

s.t. $\sum_{j=1}^n x_{ij} \leq W_i \quad \forall i$

$x_{ij} \in \mathbb{Z}_+$

where $q_{ij} = 1-p_{ij}$

The objective function becomes $\sum_{j=1}^n \prod_{i=1}^m e^{ln(q_{ij})x_{ij}}$

and then $\sum_{j=1}^n e^{\sum_{i=1}^m ln(q_{ij})x_{ij}}$.

Since $x$ is the variable, we have the sum of exponentials of a linear function of $x$ -- it's a mixed integer DCP.

Here's a quick script for it in CVXPY:

import cvxpy as cp
import numpy as np

# Define the problem data
n = 5 # number of targets
m = 3 # number of weapon types
np.random.seed(1)
q = np.random.rand(n,m)*.9 + .1 # Survival probability
W = np.random.randint(1,10,m) # Number of weapons of each type

# Define the CVXPY problem.
x = cp.Variable((n,m), integer=True)
weighted_weapons = cp.multiply(x, np.log(q)) # (n,m)
survival_probs = cp.exp(cp.sum(weighted_weapons, axis=1)) # (n,)
objective = cp.Minimize(cp.sum(survival_probs))
cons = [cp.sum(x, axis=0) <= W, x >= 0]

# Solve
prob = cp.Problem(objective, cons)
prob.solve(verbose=True)
prob.status # Optimal
prob.value # 0.0785
np.set_printoptions(precision=3, suppress=True)
x.value
# array([[0., 0., 2.],
#        [0., 3., 0.],
#        [3., 0., 0.],
#        [1., 4., 0.],
#        [4., 0., 0.]])
```