Advanced Mixed-Integer Programming (MIP) Example Using CVXPY

September 20, 2024

Advanced Mixed-Integer Programming (MIP) Example Using CVXPY

Problem: Facility Layout Problem with Sequencing and Capacity Constraints

The Facility Layout Problem is a complex optimization problem where multiple facilities need to be arranged in a layout to minimize operational costs, such as transportation between facilities, while respecting various constraints like sequencing, capacity, and adjacency requirements.

This problem becomes more challenging when sequencing and capacity constraints are introduced, turning it into a non-trivial MIP problem.

Problem Setup

Objective:

The goal is to place facilities on a grid layout to minimize the total transportation cost, considering:

Sequencing requirements (some facilities must be placed in a specific order).
Capacity constraints (some facilities can only be placed in specific regions due to space limitations).

Variables:

$(x_{ij})$: Binary variable that equals $1$ if facility $(i)$ is placed in position $(j)$, $0$ otherwise.
$(y_{ij})$: Binary variable that equals $1$ if facility $(i)$ precedes facility $(j)$ in the sequence.

Parameters:

$(c_{ij})$: Cost of placing facility $(i)$ next to facility $(j)$.
$(p_i)$: Capacity requirement for each facility.
$(g_j)$: Capacity available at each grid position.

Constraints:

Unique Placement: Each facility must be placed in exactly one position.
No Overlap: Each position can hold at most one facility.
Sequencing: Facilities must follow a predefined sequence if required.
Capacity: Facilities can only be placed in positions that meet their capacity requirements.

Mathematical Formulation

Objective:
$$
\text{Minimize} \quad \sum_{i} \sum_{j} c_{ij} x_{ij}
$$
Constraints:
1. Each facility is placed in exactly one position:
  $$
  \sum_{j} x_{ij} = 1, \quad \forall i
  $$
2. Each position can hold at most one facility:
  $$
  \sum_{i} x_{ij} \leq 1, \quad \forall j
  $$
3. Sequencing constraints between facilities:
  $$
  y_{ij} + y_{ji} = 1, \quad \forall i \neq j
  $$
4. Capacity constraints:
  $$
  p_i \times x_{ij} \leq g_j, \quad \forall i, j
  $$

Solving with CVXPY

Here’s a detailed implementation using $CVXPY$:

import cvxpy as cp
import numpy as np

# Number of facilities and grid positions
num_facilities = 4
num_positions = 4

# Cost of placing each facility next to each other (symmetric matrix)
c = np.array([
    [0, 2, 3, 4],
    [2, 0, 5, 1],
    [3, 5, 0, 6],
    [4, 1, 6, 0]
])

# Capacity requirements for each facility
p = np.array([2, 1, 3, 2])

# Available capacity at each grid position
g = np.array([3, 2, 2, 4])

# Decision variables
x = cp.Variable((num_facilities, num_positions), boolean=True)  # Placement variable
y = cp.Variable((num_facilities, num_facilities), boolean=True)  # Sequencing variable

# Objective: Minimize the total cost of placing facilities
# The objective is rewritten to respect DCP rules
objective = cp.Minimize(cp.sum(c * cp.abs(x - x.T)))

# Constraints
constraints = []

# Each facility must be placed in exactly one position
for i in range(num_facilities):
    constraints.append(cp.sum(x[i, :]) == 1)

# Each position can hold at most one facility
for j in range(num_positions):
    constraints.append(cp.sum(x[:, j]) <= 1)

# Sequencing constraints: ensure order between facilities
for i in range(num_facilities):
    for j in range(num_facilities):
        if i != j:
            constraints.append(y[i, j] + y[j, i] == 1)  # Only one direction allowed

# Capacity constraints: facility can only be placed in positions meeting capacity requirements
for i in range(num_facilities):
    for j in range(num_positions):
        constraints.append(p[i] * x[i, j] <= g[j])

# Solve the problem
problem = cp.Problem(objective, constraints)
problem.solve(solver=cp.GLPK_MI)  # Using GLPK_MI as the solver

# Results
print("Status:", problem.status)
print("Optimal Value (Total Cost):", problem.value)
print("Placement of Facilities (x):")
print(x.value)
print("Sequencing of Facilities (y):")
print(y.value)

Explanation of the Code:

Decision Variables:
- x[i, j]: Binary variable indicating whether facility $( i )$ is placed in position $( j )$.
- y[i, j]: Binary variable indicating whether facility $( i )$ precedes facility $( j )$ in the sequence.
Objective:
- The objective is to minimize the total placement cost, computed by multiplying the placement matrix with the cost matrix.
Constraints:
- Each facility must be assigned exactly one position, and each position can hold at most one facility.
- Sequencing constraints ensure that the required order between facilities is maintained.
- Capacity constraints enforce that a facility can only be placed in a position with sufficient available capacity.
Solver:
- The problem is solved using GLPK_MI, an open-source solver suitable for MIP problems.

Output:

Status: optimal
Optimal Value (Total Cost): 0.0
Placement of Facilities (x):
[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 0. 1.]
 [0. 0. 1. 0.]]
Sequencing of Facilities (y):
[[0. 1. 1. 1.]
 [0. 0. 1. 1.]
 [0. 0. 0. 1.]
 [0. 0. 0. 0.]]

Conclusion:

This advanced Mixed-Integer Programming (MIP) example demonstrates how to solve a complex facility layout problem with sequencing and capacity constraints using $CVXPY$.

The use of binary variables allows us to model the discrete nature of facility placements and sequencing requirements, making MIP a powerful tool for such complex layout optimization problems.

Second-Order Cone Programming (SOCP) Example Using CVXPY

September 19, 2024

Second-Order Cone Programming (SOCP) Example Using CVXPY

Problem: Portfolio Optimization with Transaction Costs

In portfolio optimization, we aim to allocate a given budget across various assets to maximize returns while controlling for risk.

In this example, we will add transaction costs (e.g., costs incurred when buying and selling assets) to the optimization problem.
This problem is naturally modeled as a Second-Order Cone Programming (SOCP) problem.

Problem Setup:

Variables:
- $( x_i )$: The amount invested in asset $( i )$.
- $( r_i )$: The expected return for asset $( i )$.
- $( \Sigma )$: The covariance matrix representing the risk between the assets.
Objective:
Maximize the expected return while controlling for risk (measured as portfolio variance) and including transaction costs.
Formulation:
$$
\text{Maximize} \quad r^T x - \gamma \cdot \sqrt{x^T \Sigma x} - \lambda \cdot \sum_i |x_i - x_i^0|
$$
where:
- $( r^T x )$ is the expected return,
- $( \gamma )$ is the risk aversion parameter,
- $( \sqrt{x^T \Sigma x} )$ is the portfolio variance (risk),
- $( \lambda )$ is the transaction cost parameter,
- $( |x_i - x_i^0| )$ is the absolute change in position (transaction cost), where $( x_i^0 )$ represents the initial positions.
Constraints:
1. Total investment must equal the available budget.
2. The portfolio must stay within certain bounds for each asset.

SOCP Formulation:

We can represent this optimization problem using SOCP constraints by rewriting the risk term $( \sqrt{x^T \Sigma x} )$ as a second-order cone constraint.

Code Implementation Using CVXPY

import cvxpy as cp
import numpy as np

# Data for the problem
n = 5  # Number of assets
np.random.seed(0)

# Expected returns for each asset
r = np.random.rand(n)

# Covariance matrix (symmetric positive semidefinite)
Sigma = np.random.randn(n, n)
Sigma = Sigma.T @ Sigma  # Make the covariance matrix positive semi-definite

# Initial positions (already held assets)
x0 = np.random.rand(n)

# Parameters
gamma = 0.1  # Risk aversion parameter
lambda_ = 0.05  # Transaction cost parameter
budget = 1  # Total budget

# Decision variable: investment in each asset
x = cp.Variable(n)

# Objective: maximize expected return - risk - transaction costs
# Risk term: sqrt(x.T * Sigma * x) is modeled as second-order cone
risk = cp.quad_form(x, Sigma)  # Quadratic form for risk (x.T @ Sigma @ x)
transaction_costs = cp.norm1(x - x0)  # l1 norm for transaction costs

# Maximize expected return minus risk and transaction costs
objective = cp.Maximize(r.T @ x - gamma * cp.norm(risk, 2) - lambda_ * transaction_costs)

# Constraints
constraints = [
    cp.sum(x) == budget,  # Total investment must equal budget
    x >= 0,  # No short-selling (no negative positions)
    x <= 1   # No asset can take more than 100% of the budget
]

# Define and solve the problem
problem = cp.Problem(objective, constraints)
problem.solve()

# Results
print("Status:", problem.status)
print("Optimal Portfolio:", x.value)
print("Optimal Value:", problem.value)

Explanation of the Code:

Data:
- r: A vector of expected returns for each asset.
- Sigma: The covariance matrix representing the risk between assets.
- x0: The initial asset positions before the optimization.
Objective:
- The objective is to maximize the expected return, minus the risk (modeled as a second-order cone) and transaction costs (modeled as the $( l_1 )$-norm of the change in positions).
- cp.quad_form(x, Sigma) represents the quadratic risk term, $( x^T \Sigma x )$.
- The $( l_1 )$-norm of x - x0 captures the absolute transaction costs.
Constraints:
- The sum of the asset allocations must equal the available budget.
- We enforce non-negative asset allocations (no short-selling) and limit the allocation to a maximum of $100$% per asset.
SOCP Structure:
- The risk term involves a second-order cone constraint because of the quadratic term $( \sqrt{x^T \Sigma x} )$, which $CVXPY$ automatically handles as an SOCP.

Output Example:

Status: optimal
Optimal Portfolio: [4.39530365e-01 3.32443990e-01 5.68332295e-10 2.28025644e-01
 1.65363907e-10]
Optimal Value: 0.49075615889954766

In this output:

The Optimal Portfolio gives the asset allocation that maximizes the return while minimizing risk and transaction costs.
The Optimal Value represents the objective value at the solution.

Conclusion:

This example demonstrates how to solve a portfolio optimization problem with transaction costs using Second-Order Cone Programming (SOCP) in $CVXPY$.

SOCP allows us to model the risk term (which involves a quadratic form) as a second-order cone, and $CVXPY$ simplifies the solution process.

This approach is particularly useful for finance problems where risk and transaction costs need to be balanced.

Mixed-Integer Programming (MIP) Example Using CVXPY

September 18, 2024

Mixed-Integer Programming (MIP) Example Using CVXPY

Problem: Facility Location Problem

The facility location problem is a classic optimization problem where a company must decide where to open facilities (warehouses) to serve customers.

The goal is to minimize the total cost, which includes fixed costs for opening facilities and transportation costs for serving customers.

Problem Setup:

Variables:
- $( x_{ij} )$: Binary decision variable, $1$ if customer $( j )$ is served by facility $( i )$, $0$ otherwise.
- $( y_i )$: Binary decision variable, $1$ if facility $( i )$ is open, $0$ otherwise.
Parameters:
- $( f_i )$: Fixed cost of opening facility $( i )$.
- $( c_{ij} )$: Transportation cost for serving customer $( j )$ from facility $( i )$.
Objective:
Minimize the total cost:
$$
\text{Minimize} \quad \sum_{i} f_i y_i + \sum_{i}\sum_{j} c_{ij} x_{ij}
$$
Constraints:
1. Each customer must be assigned to exactly one facility:
  $$
  \sum_{i} x_{ij} = 1, \quad \forall j
  $$
2. A customer can only be served by an open facility:
  $$
  x_{ij} \leq y_i, \quad \forall i, j
  $$
3. Binary constraints on the variables:
  $$
  x_{ij} \in {0, 1}, \quad y_i \in {0, 1}, \quad \forall i, j
  $$

Solving with CVXPY:

We can solve this problem using $CVXPY$, a Python library for convex optimization, which also supports Mixed-Integer Programming.

Here’s a step-by-step guide with Python code.

import cvxpy as cp
import numpy as np

# Data (example)
num_facilities = 3  # Number of facilities
num_customers = 5   # Number of customers

# Fixed costs for opening each facility
f = np.array([100, 200, 300])

# Transportation costs from each facility to each customer
c = np.array([
    [8, 10, 7, 6, 9],
    [6, 8, 5, 7, 10],
    [10, 7, 6, 8, 7]
])

# Decision variables
# x_ij: binary, 1 if customer j is served by facility i
x = cp.Variable((num_facilities, num_customers), boolean=True)

# y_i: binary, 1 if facility i is open
y = cp.Variable(num_facilities, boolean=True)

# Objective: Minimize fixed and transportation costs
# Updated to use elementwise multiplication (cp.multiply)
objective = cp.Minimize(cp.sum(f @ y) + cp.sum(cp.multiply(c, x)))

# Constraints
constraints = []

# Each customer must be assigned to exactly one facility
for j in range(num_customers):
    constraints.append(cp.sum(x[:, j]) == 1)

# Customers can only be assigned to open facilities
for i in range(num_facilities):
    for j in range(num_customers):
        constraints.append(x[i, j] <= y[i])

# Solve the problem using another solver such as GLPK_MI
problem = cp.Problem(objective, constraints)
problem.solve(solver=cp.GLPK_MI)

# Results
print("Status:", problem.status)
print("Total Cost:", problem.value)
print("Facilities to Open:", y.value)
print("Customer Assignments:")
print(x.value)

Explanation:

Decision Variables:
- x: This is a ( num_facilities $\times$ num_customers ) matrix where each entry represents whether a customer is served by a facility.
  We use boolean=True to enforce the binary constraint.
- y: This is a vector of length num_facilities that indicates whether a facility is open.
Objective:
- The objective is to minimize the sum of the fixed costs $( f_i )$ for opening facilities and the transportation costs $( c_{ij} )$ for serving customers.
Constraints:
- We enforce the constraint that each customer is assigned to exactly one facility using cp.sum(x[:, j]) == 1.
- We ensure that customers are only assigned to open facilities with x[i, j] <= y[i].
Solving:
- We use the GLPK_MI solver (solver=cp.GLPK_MI), which is well-suited for mixed-integer programming problems.
- After solving, we can inspect the status, total cost, which facilities to open, and the customer assignments.

Output:

Status: optimal
Total Cost: 140.0
Facilities to Open: [1. 0. 0.]
Customer Assignments:
[[1. 1. 1. 1. 1.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]

This indicates that:

Facilities $1$ should be opened.
Customers $1$, $2$, $3$, $4$, and $5$ are served by facility 1.

Conclusion:

This example demonstrates how to model and solve a Mixed-Integer Programming (MIP) problem using $CVXPY$.

The facility location problem is a common application of MIP, and it shows how we can optimize costs by using binary decision variables to represent discrete choices (e.g., opening facilities).

Robust Optimization Example Using CVXPY

September 17, 2024

Robust Optimization Example Using CVXPY

Robust optimization is a technique used to handle uncertainty in optimization problems.

In contrast to standard optimization, where we assume exact data, robust optimization considers that some of the data (like costs, constraints, or coefficients) are uncertain and vary within known bounds.

This ensures that the solution is feasible and performs well even in the worst-case scenario of these uncertainties.

In this example, we will solve a robust portfolio optimization problem where we aim to minimize the worst-case risk (variance) of a portfolio under uncertain returns.

Problem Description: Robust Portfolio Optimization

Imagine we have an investor who wants to allocate their capital across several assets.

The investor is uncertain about the expected returns of these assets but knows that the returns lie within a given uncertainty range.
The goal is to find an optimal portfolio allocation that minimizes the worst-case risk (variance) while ensuring the expected return exceeds a certain threshold.

Variables:

$( x )$: Vector of asset weights (the proportion of capital allocated to each asset).

Objective:

Minimize the worst-case variance of the portfolio.

Constraints:

The total portfolio allocation should sum to $1$ (i.e., $( \sum x_i = 1 )$).
The expected return should be at least a target value $( R_{\text{min}} )$.
The asset returns are uncertain but lie within a given uncertainty set.

We will now define and solve this robust optimization problem in $Python$ using $CVXPY$.

Step-by-Step Solution

import cvxpy as cp
import numpy as np

# Step 1: Define the problem data
n = 4  # Number of assets
np.random.seed(42)

# Covariance matrix of returns (representing risk/variance)
Sigma = np.array([
    [0.1, 0.02, 0.04, 0.01],
    [0.02, 0.08, 0.03, 0.02],
    [0.04, 0.03, 0.09, 0.03],
    [0.01, 0.02, 0.03, 0.07]
])

# Expected returns (uncertain)
mu = np.array([0.1, 0.12, 0.14, 0.08])  # Nominal expected returns
delta_mu = np.array([0.02, 0.03, 0.02, 0.01])  # Uncertainty in returns

# Minimum target return
R_min = 0.11

# Step 2: Define the decision variables
x = cp.Variable(n)  # Portfolio weights

# Step 3: Define the robust constraints and objective

# Constraint 1: The total portfolio allocation should sum to 1
constraints = [cp.sum(x) == 1, x >= 0]

# Robust constraint: Expected return should be at least R_min for the worst-case returns
worst_case_return = mu - delta_mu  # Worst-case scenario for returns
constraints.append(worst_case_return @ x >= R_min)

# Objective: Minimize the worst-case portfolio variance
portfolio_variance = cp.quad_form(x, Sigma)
objective = cp.Minimize(portfolio_variance)

# Step 4: Define and solve the optimization problem
problem = cp.Problem(objective, constraints)
problem.solve()

# Step 5: Output the results
print("Status:", problem.status)
print("Optimal portfolio allocation:", x.value)
print("Minimum worst-case risk (variance):", portfolio_variance.value)

Explanation

Problem Data:
- We define a covariance matrix Sigma, representing the variance and correlation between the returns of the assets. The diagonal elements represent the variance of individual assets, while the off-diagonal elements represent the covariance between pairs of assets.
- The vector mu contains the nominal (expected) returns for each asset, and delta_mu represents the uncertainty in these returns.
- The target return $( R_{\text{min}} )$ is set to $0.11$, meaning the investor wants the portfolio to yield at least an $11$% return, even in the worst-case scenario.
Decision Variables:
- The variable x represents the portfolio weights, i.e., the proportion of the total capital allocated to each asset.
Constraints:
- The total allocation must sum to $1$ (cp.sum(x) == 1), meaning all the capital must be invested.
- The worst-case expected return (worst_case_return @ x) should be at least $( R_{\text{min}} )$, ensuring the portfolio remains profitable under uncertainty.
- We also require non-negative portfolio weights (x >= 0), assuming no short selling.
Objective:
- The objective is to minimize the portfolio’s worst-case variance (cp.quad_form(x, Sigma)), which is a quadratic form representing the portfolio’s risk.
Solution:
- We solve the problem using problem.solve(). If the problem is feasible, it will return the optimal portfolio weights that minimize the worst-case risk while satisfying the constraints.

Output

1
2
3

Status: optimal
Optimal portfolio allocation: [0.02536232 0.29347826 0.67753623 0.00362319]
Minimum worst-case risk (variance): 0.062065217391304356

Interpretation of Results

Optimal Portfolio Allocation:
The optimal portfolio weights allocate the investor’s capital across the four assets.
For instance, $25$% of the total capital is allocated to the first asset, $29$% to the second, and so on.
Minimum Worst-Case Risk:
The worst-case risk (variance) of the portfolio is minimized at around $0.062$, indicating the portfolio has been optimized to minimize risk under uncertain returns.

Conclusion

This example demonstrated how to use $CVXPY$ to solve a robust optimization problem in the context of portfolio optimization.

Robust optimization allows us to account for uncertainty in the data (such as uncertain returns in this case) and find solutions that are robust to worst-case scenarios.

The problem ensures that the portfolio achieves the target return even in the worst case while minimizing the overall risk.

This technique is useful in finance and other fields where uncertainty is present, and decision-makers seek solutions that perform well under a range of possible conditions.

Least Squares and Regression with CVXPY

September 16, 2024

Least Squares and Regression with CVXPY

Least squares is a popular method for solving regression problems, where we want to fit a model to data by minimizing the sum of the squared differences (errors) between the observed values and the values predicted by the model.

In this example, we will demonstrate how to use $CVXPY$ to solve a $linear$ $regression$ problem using least squares.

Problem Description: Linear Regression

We are given a dataset consisting of several observations.
Each observation includes:

A set of features (independent variables).
A target value (dependent variable).

Our goal is to find a linear relationship between the features and the target, meaning we want to find the best-fitting line (or hyperplane) that predicts the target based on the features.

Mathematically, we can express this as:
$$
y = X\beta + \epsilon
$$
Where:

$( X )$ is a matrix of features (each row corresponds to an observation and each column to a feature).
$( y )$ is a vector of observed target values.
$( \beta )$ is the vector of unknown coefficients we want to estimate.
$( \epsilon )$ represents the errors (differences between the predicted and observed values).

We want to find $( \beta )$ that minimizes the sum of squared errors:
$$
\min_\beta | X\beta - y |_2^2
$$

Step-by-Step Solution with CVXPY

Here is how to solve this least squares problem using $CVXPY$:

import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt

# Step 1: Generate synthetic data for a regression problem
np.random.seed(1)
n = 50  # Number of data points
m = 3   # Number of features

# Randomly generate feature matrix X and true coefficients beta_true
X = np.random.randn(n, m)
beta_true = np.array([3, -1, 2])  # True coefficients
y = X @ beta_true + 0.5 * np.random.randn(n)  # Observed values with some noise

# Step 2: Define the variable (beta) for the regression coefficients
beta = cp.Variable(m)

# Step 3: Define the objective function (least squares error)
objective = cp.Minimize(cp.sum_squares(X @ beta - y))

# Step 4: Define and solve the problem
problem = cp.Problem(objective)
problem.solve()

# Step 5: Output the results
print("Status:", problem.status)
print("Estimated coefficients (beta):", beta.value)
print("True coefficients:", beta_true)

# Step 6: Visualize the data and the fitted values
y_pred = X @ beta.value
plt.scatter(range(n), y, color='blue', label='Observed')
plt.scatter(range(n), y_pred, color='red', label='Predicted')
plt.xlabel('Data Point')
plt.ylabel('Target Value')
plt.title('Least Squares Regression Fit')
plt.legend()
plt.show()

Detailed Explanation:

Data Generation:
In this example, we generate synthetic data for simplicity.
The matrix $( X )$ contains $50$ data points with $3$ features, and we generate the target values $( y )$ using a known set of coefficients $( \beta_{\text{true}} = [3, -1, 2] )$ plus some noise to simulate real-world observations.
Decision Variable:
The unknowns are the regression coefficients $( \beta )$, represented as a $CVXPY$ variable beta of size $3$ (since we have $3$ features).
Objective Function:
The goal is to minimize the sum of squared errors between the observed target values $( y )$ and the predicted values $( X\beta )$.
In $CVXPY$, this is represented by the function cp.sum_squares(X @ beta - y), which computes the sum of the squared residuals.
Problem Definition and Solution:
The cp.Problem function defines the least squares optimization problem, and problem.solve() finds the optimal solution for $( \beta )$.
Results:
After solving the problem, the estimated coefficients $( \beta )$ are printed.
These should be close to the true coefficients used to generate the data.
We also visualize the observed values versus the predicted values to show how well the model fits the data.

Output:

Estimated Coefficients: The estimated coefficients $( \hat{\beta} )$ are close to the true coefficients used to generate the data ($( \beta_{\text{true}} = [3, -1, 2] )$), which indicates a good fit.

Visualization of Results:

The following plot shows the observed values (in blue) and the predicted values (in red).
A good fit would show the predicted values closely following the observed ones.

Blue points: Actual (observed) target values.
Red points: Predicted target values based on the fitted model.

Interpretation:

Minimizing the Sum of Squared Errors:
The least squares method minimizes the squared differences between the predicted and observed values.
This produces the “best fit” line for the given data, in the sense that the total prediction error is minimized.
Optimal Solution:
Since the status of the problem is optimal, $CVXPY$ successfully found the regression coefficients that minimize the objective function.

Conclusion:

In this example, we demonstrated how to solve a linear regression problem using least squares with $CVXPY$.
The least squares method is widely used in machine learning and statistics for fitting linear models to data.

By minimizing the sum of squared errors, we can estimate the coefficients that best explain the relationship between the features and the target values.

Quadratic Programming Example Using CVXPY

September 15, 2024

Quadratic Programming Example Using CVXPY

$Quadratic Programming (QP)$ is an optimization problem where the objective function is quadratic, and the constraints are linear.

A typical $QP$ problem can be written as:

$$
\text{Minimize } \frac{1}{2} x^T Q x + c^T x
$$
subject to:
$$
Ax \leq b \quad \text{and} \quad Ex = d
$$

Where:

$( Q )$ is a positive semi-definite matrix (for convexity),
$( c )$ is a vector of linear coefficients,
$( A )$, $( b )$, $( E )$, and $( d )$ represent the inequality and equality constraints, respectively.

Example: Portfolio Optimization

We will solve a portfolio optimization problem using $QP$.

The goal is to minimize the risk (variance) of a portfolio, subject to achieving a certain return and satisfying investment constraints.

Problem Formulation:

Objective: Minimize the portfolio’s risk, measured as the variance of returns.
Constraints:
- The total weight of assets must sum to $1$ (fully invested).
- Each weight must be non-negative (no short selling).
- The portfolio’s expected return must be at least a target return.

Data Setup:

The problem has 3 assets with the following parameters:
- Expected returns: $([0.1, 0.2, 0.15])$
- Covariance matrix (risk):
  $$
  \begin{pmatrix}
  0.005 & -0.010 & 0.004 \
  -0.010 & 0.040 & -0.002 \
  0.004 & -0.002 & 0.023 \
  \end{pmatrix}
  $$
- Target return: $(0.16)$

Step-by-Step Solution:

Define the Variables and Problem:
Let $( x = \begin{bmatrix} x_1 \ x_2 \ x_3 \end{bmatrix} )$ represent the weights of the three assets in the portfolio.
Objective:
The objective is to minimize the variance (quadratic term):
$$
\frac{1}{2} x^T Q x
$$
where $( Q )$ is the covariance matrix of the asset returns.
Constraints:
- The sum of portfolio weights must be 1: $( x_1 + x_2 + x_3 = 1 )$
- The expected return of the portfolio must be at least $(0.16)$:
$$
0.1x_1 + 0.2x_2 + 0.15x_3 \geq 0.16
$$
- No short selling: $( x_1, x_2, x_3 \geq 0 )$

Python Implementation:

We will use the $cvxpy$ library to solve this $QP$ problem.

import cvxpy as cp
import numpy as np

# Step 1: Define the problem data
Q = np.array([[0.005, -0.010, 0.004],
              [-0.010, 0.040, -0.002],
              [0.004, -0.002, 0.023]])  # Covariance matrix

expected_returns = np.array([0.1, 0.2, 0.15])  # Expected returns for the assets
target_return = 0.16  # Target portfolio return

# Step 2: Define the decision variable (portfolio weights)
x = cp.Variable(3)

# Step 3: Define the objective function (minimize variance)
objective = cp.Minimize((1/2) * cp.quad_form(x, Q))

# Step 4: Define the constraints
constraints = [
    cp.sum(x) == 1,  # Fully invested (weights sum to 1)
    expected_returns @ x >= target_return,  # Target return constraint
    x >= 0  # No short selling constraint
]

# Step 5: Formulate the problem
problem = cp.Problem(objective, constraints)

# Step 6: Solve the problem
problem.solve()

# Step 7: Output results
print("Status:", problem.status)
print("Optimal portfolio allocation:", x.value)
print("Minimum risk (variance):", problem.value)

Detailed Explanation:

Covariance Matrix ( Q ):
This matrix represents the risk (variance and covariance) associated with the assets in the portfolio.
The term cp.quad_form(x, Q) calculates the quadratic form $( x^T Q x )$, representing the portfolio variance.
Objective Function:
The objective is to minimize the portfolio variance.
In cvxpy, we use the cp.Minimize() function, where the quadratic form represents the risk.
Constraints:
- Fully Invested: cp.sum(x) == 1 ensures that the portfolio is fully invested (the weights sum to 1).
- Target Return: expected_returns @ x >= target_return ensures that the portfolio’s expected return is at least $0.16$.
- No Short Selling: x >= 0 ensures that no asset has a negative weight (i.e., no short selling).
Solve the Problem:
The problem.solve() function solves the $QP$ problem and finds the optimal portfolio allocation that minimizes risk while satisfying all the constraints.

Output:

1
2
3

Status: optimal
Optimal portfolio allocation: [0.26055046 0.46055046 0.27889908]
Minimum risk (variance): 0.004140183486238534

Interpretation of Results:

Optimal Portfolio Allocation:
The optimal weights for the assets are approximately $26.06$% in Asset $1$, $46.06$% in Asset $2$, and $27.89$% in Asset $3$.
These weights minimize the portfolio variance while achieving the target return.
Minimum Risk:
The minimum variance (risk) is approximately $0.00414$.
This represents the lowest possible portfolio risk that satisfies the return constraint.

Conclusion:

This example demonstrates how to solve a quadratic programming problem using Python’s $cvxpy$ library.

The problem was formulated as a portfolio optimization task, where we minimized risk subject to return and no-short-selling constraints.

Maximizing the Minimum Eigenvalue of a Matrix Using Semidefinite Programming in CVXPY

September 14, 2024

Maximizing the Minimum Eigenvalue of a Matrix Using Semidefinite Programming in CVXPY

Semidefinite Programming (SDP) in CVXPY

$Semidefinite programming (SDP)$ is an optimization problem where the decision variable is a matrix, and the constraints include that this matrix must be positive semidefinite.

$SDP$ problems arise in various fields, such as control theory, signal processing, quantum mechanics, and machine learning.

Problem Setup:

The general form of a semidefinite program is:

$$
\text{minimize } \quad \text{tr}(CX)
$$
$$
\text{subject to } \quad \text{tr}(A_iX) = b_i \quad \text{for each } i
$$
$$
X \succeq 0
$$

Where:

$(X)$ is the decision variable (a symmetric matrix).
$(\text{tr}(C X))$ is the objective function (trace of matrix multiplication).
The matrix $(X \succeq 0)$ indicates that $(X)$ is positive semidefinite, meaning all its eigenvalues are non-negative.

Example: Maximizing the Minimum Eigenvalue of a Matrix

Let’s solve an $SDP$ example where we want to maximize the minimum eigenvalue of a symmetric matrix $( X )$, subject to some constraints.

Problem:

Given a symmetric matrix $( X )$, we want to:

Maximize the minimum eigenvalue of $( X )$.
Ensure that $( X )$ satisfies some linear constraints.

Formulation:

We want to maximize $( \lambda_{\text{min}}(X) )$, the smallest eigenvalue of $( X )$, which can be written as an $SDP$ using a slack variable $( t )$ and ensuring $( X - tI \succeq 0 )$ (i.e., the matrix $( X - tI )$ is positive semidefinite).

Step-by-Step Solution using CVXPY:

import cvxpy as cp
import numpy as np

# Step 1: Define the size of the matrix
n = 3  # Size of the matrix X (3x3 matrix in this example)

# Step 2: Define the decision variable (X is a symmetric matrix)
X = cp.Variable((n, n), symmetric=True)

# Step 3: Define the slack variable (t) representing the minimum eigenvalue
t = cp.Variable()

# Step 4: Define the objective function (maximize t)
# This corresponds to maximizing the smallest eigenvalue
objective = cp.Maximize(t)

# Step 5: Define the constraints
# The first constraint is X - tI >= 0 (positive semidefinite)
constraints = [X - t * np.eye(n) >> 0]

# Add some additional constraints to make the problem non-trivial
# For example, trace(X) = 10, and X[0, 1] = 1 (arbitrary example constraints)
constraints += [cp.trace(X) == 10]
constraints += [X[0, 1] == 1]

# Step 6: Formulate the problem
problem = cp.Problem(objective, constraints)

# Step 7: Solve the problem
problem.solve()

# Step 8: Output results
print("Status:", problem.status)
print("Optimal value of t (minimum eigenvalue):", t.value)
print("Optimal solution for X:\n", X.value)

Explanation:

Step-by-Step Breakdown:

Decision Variable:
- We define X as a symmetric matrix variable of size $(n \times n)$.
  In our case, we chose $(n = 3)$, so $(X)$ is a $(3 \times 3)$ symmetric matrix.
Slack Variable $(t)$:
- The variable t represents the smallest eigenvalue of the matrix $(X)$.
  The objective is to maximize $(t)$, which corresponds to maximizing the minimum eigenvalue of $(X)$.
Objective Function:
- The goal is to maximize $(t)$, i.e., maximize the minimum eigenvalue of $(X)$.
Constraints:
- Semidefinite Constraint: $(X - tI \succeq 0)$ (i.e., the matrix $(X - tI)$ must be positive semidefinite). This ensures that $(t)$ is a lower bound on the eigenvalues of $(X)$.
  All eigenvalues of $(X)$ must be greater than or equal to $(t)$.
- Additional Constraints:
  - We impose $( \text{tr}(X) = 10 )$ (i.e., the sum of the diagonal elements of $(X) equals 10)$.
    This is an arbitrary constraint to ensure the matrix is not trivial.
  - We also constrain $( X[0, 1] = 1 )$, just to add more structure to the problem.
Problem Definition:
- The problem is defined using cp.Problem(), with the objective being the maximization of$ (t)$ and the constraints defined above.
Solve:
- The solve() method is called to solve the semidefinite program.
Results:
- The status of the solution is printed (Optimal, Infeasible, or Unbounded).
- The optimal value of $(t)$, which corresponds to the maximum minimum eigenvalue of (X), is printed.
- The optimal solution matrix $(X)$ is also displayed.

Output:

Status: optimal
Optimal value of t (minimum eigenvalue): 2.6666580666942266
Optimal solution for X:
 [[3.666671   1.00000001 0.        ]
 [1.00000001 3.666671   0.        ]
 [0.         0.         2.666658  ]]

Optimal value of $(t)$: The minimum eigenvalue of the matrix $(X)$ is approximately $2.67$.
Optimal solution for $(X)$: The matrix $(X)$ that satisfies the given constraints and maximizes the minimum eigenvalue is shown.

Explanation of Results:

Semidefinite Constraint: The constraint $$( X - tI \succeq 0 )$$ ensures that all eigenvalues of $(X)$ are greater than or equal to $(t)$.
Thus, maximizing $(t)$ corresponds to maximizing the smallest eigenvalue of $(X)$.
Trace Constraint: The sum of the diagonal elements of $(X)$ is constrained to be 10.
Off-Diagonal Constraint: The constraint $( X[0, 1] = 1 )$ ensures that there’s a relationship between the elements of $(X)$, adding complexity to the optimization.

Use Cases of Semidefinite Programming:

Control Systems: $SDP$ can be used in control theory to design controllers that guarantee system stability by finding matrices that satisfy Lyapunov stability criteria.
Structural Design: Optimizing stiffness matrices for structures to maximize strength while minimizing material usage.
Quantum Information: $SDP$ is used to analyze quantum states, where matrices representing states must be positive semidefinite.
Finance: Portfolio optimization and risk management, where covariance matrices must satisfy certain positive semidefinite constraints.

Conclusion:

This example demonstrates how $CVXPY$ can be used to solve semidefinite programming problems by formulating the problem with matrix variables, semidefinite constraints, and other linear constraints.

The power of $SDP$ lies in its ability to optimize over matrices while ensuring positive semidefiniteness, making it a crucial tool in many advanced applications across different fields.

Minimizing Energy Consumption in Mass Transport Using PuLP Optimization

September 13, 2024

Minimizing Energy Consumption in Mass Transport Using PuLP Optimization

Here’s an example of using $PuLP$ to solve a physics-based optimization problem.

Let’s solve the “minimum energy transportation problem”, where we aim to transport a mass from one point to another, minimizing the total energy consumed.

Problem Statement:

We want to move a mass between several locations while minimizing the energy spent.

The energy required to move the mass depends on the distance between locations and the work done against gravity (if there are changes in elevation).

Given:

A set of locations with known distances and elevation changes.
A mass to transport.
We need to determine the optimal path to minimize the energy required.

Assumptions:

The total energy required consists of two components: kinetic energy (from horizontal distance) and potential energy (from elevation changes).
We’ll use simple physics formulas:
- Potential Energy: $( E_p = mgh ) $
- Kinetic Energy: $( E_k = \frac{1}{2} m v^2 )$, but we will simplify this to a linear relation for distance.

Solving the Problem with PuLP:

import pulp

# Step 1: Define problem data
locations = ['A', 'B', 'C', 'D']
distances = {('A', 'B'): 10, ('A', 'C'): 20, ('B', 'C'): 15, ('B', 'D'): 30, ('C', 'D'): 25}
elevations = {'A': 0, 'B': 100, 'C': 50, 'D': 200}  # in meters
mass = 10  # mass of object (in kg)
gravity = 9.8  # gravitational constant (m/s^2)

# Step 2: Define the optimization problem
prob = pulp.LpProblem("Minimize_Energy", pulp.LpMinimize)

# Step 3: Define decision variables (binary: 1 if we move from location i to j, 0 otherwise)
x = pulp.LpVariable.dicts("path", distances, 0, 1, pulp.LpBinary)

# Step 4: Define the objective function (minimize total energy)
energy_cost = pulp.lpSum(
    [x[i, j] * (mass * gravity * abs(elevations[j] - elevations[i]) + distances[i, j]) for (i, j) in distances]
)
prob += energy_cost, "Total_Energy"

# Step 5: Constraints
# Start at A and end at D
prob += pulp.lpSum([x['A', j] for j in locations if ('A', j) in distances]) == 1, "Start_at_A"
prob += pulp.lpSum([x[i, 'D'] for i in locations if (i, 'D') in distances]) == 1, "End_at_D"

# Each intermediate location must be left if it's visited
for loc in locations:
    if loc not in ['A', 'D']:
        prob += pulp.lpSum([x[i, loc] for i in locations if (i, loc) in distances]) == \
                pulp.lpSum([x[loc, j] for j in locations if (loc, j) in distances]), f"flow_conservation_{loc}"

# Step 6: Solve the problem
prob.solve()

# Step 7: Output results
print("Status:", pulp.LpStatus[prob.status])

print("Optimal Path:")
for i, j in distances:
    if pulp.value(x[i, j]) == 1:
        print(f"Move from {i} to {j}")

print("Minimum Energy Required:", pulp.value(prob.objective))

Explanation:

Locations and Distances: We define a set of locations $(A, B, C, D)$ and the distances between them.
We also define the elevation of each location.
Energy Calculation: The total energy required consists of two parts:
- The potential energy due to elevation changes $(E_p = mgh)$.
- The kinetic energy due to horizontal movement is simplified by adding distance as a cost.
Objective Function: We minimize the total energy required to transport the object from the start to the end location.
Constraints:
- The mass starts at location $A$ and must end at location $D$.
- There are flow conservation constraints to ensure that if the mass enters an intermediate location, it also leaves it.
Solution: $PuLP$ solves the problem, providing the optimal path that minimizes energy consumption.

Output:

Status: Optimal
Optimal Path:
Move from A to B
Move from B to D
Minimum Energy Required: 19640.0

Use Case:

This kind of problem can be used in robotics, logistics, or space missions, where minimizing energy usage is crucial for efficient operation, especially in systems dealing with gravity and elevation changes.

Production Planning Example with PuLP

September 12, 2024

Production Planning Example with PuLP

Let’s consider a simple Production Planning problem for a factory that produces two products: Product A and Product B.

The factory has limited resources: labor hours and raw material.

The goal is to maximize the profit while respecting the constraints of available resources.

Problem Statement:

Product A requires:
- $2$ labor hours per unit
- $3$ kg of raw material per unit
- Yields a profit of $30 per unit
Product B requires:
- $3$ labor hours per unit
- $1$ kg of raw material per unit
- Yields a profit of $20 per unit

The factory has:

$100$ labor hours available per week
$90$ kg of raw material available per week

The factory needs to decide how many units of each product to produce to maximize profit, given the constraints.

Mathematical Model:

Let:

$(x_1)$ = number of units of Product A produced
$(x_2)$ = number of units of Product B produced

Objective:
$$
\text{Maximize } 30x_1 + 20x_2
$$
Subject to:
$$
2x_1 + 3x_2 \leq 100 \quad \text{(labor hours constraint)}
$$
$$
3x_1 + 1x_2 \leq 90 \quad \text{(raw material constraint)}
$$
$$
x_1 \geq 0, \quad x_2 \geq 0 \quad \text{(non-negativity constraint)}
$$

Solution using PuLP:

import pulp

# Define the problem
problem = pulp.LpProblem("Production_Planning", pulp.LpMaximize)

# Decision variables
x1 = pulp.LpVariable('x1', lowBound=0, cat='Continuous')  # Units of Product A
x2 = pulp.LpVariable('x2', lowBound=0, cat='Continuous')  # Units of Product B

# Objective function
problem += 30 * x1 + 20 * x2, "Profit"

# Constraints
problem += 2 * x1 + 3 * x2 <= 100, "Labor_hours"
problem += 3 * x1 + 1 * x2 <= 90, "Raw_material"

# Solve the problem
problem.solve()

# Output the results
print("Status:", pulp.LpStatus[problem.status])
print(f"Optimal number of Product A to produce: {x1.varValue}")
print(f"Optimal number of Product B to produce: {x2.varValue}")
print(f"Maximum Profit: ${pulp.value(problem.objective)}")

Explanation:

Decision variables represent the number of units of each product to be produced.
Objective function aims to maximize the total profit from producing both products.
Constraints ensure that the production doesn’t exceed the available labor hours and raw material.

Results:

Status: Optimal
Optimal number of Product A to produce: 24.285714
Optimal number of Product B to produce: 17.142857
Maximum Profit: $1071.4285599999998

Explanation of the Results:

The problem has been solved, and the solution status is “$Optimal$”, meaning the solution found satisfies all the constraints and maximizes the profit under the given conditions.

Optimal number of Product A to produce: $24.29$ units (rounded to two decimal places)
Optimal number of Product B to produce: $17.14$ units (rounded to two decimal places)

Since the solution involves fractional units (which might represent part of a batch or a continuous process in certain industries), the optimal production mix is not whole numbers.

Maximum Profit: $1071.43

This means that producing approximately $24.29$ units of Product A and $17.14$ units of Product B will yield the highest possible profit of $1071.43, while respecting the constraints on labor hours and raw materials.

If the production process must involve whole units, you might consider rounding these values or solving the problem using integer programming, depending on the context.

Solving the Diet Problem with PuLP

September 11, 2024

An Example of Cost-Effective Nutrition Optimization

Example: Diet Problem

The diet problem is a classic $optimization$ $problem$ that aims to find the most $cost$-$effective$ way to meet nutritional requirements using a selection of foods.

Let’s solve a simple example using the $PuLP$ library in $Python$.

Problem Statement:

You are trying to meet your daily nutritional needs using three different foods.

Each food has a $cost$ and provides a certain amount of $calories$, $protein$, and $fat$.

Your goal is to minimize the cost while ensuring that you meet the minimum daily requirements for $calories$, $protein$, and $fat$.

Data:

Foods: Chicken, Beef, Vegetables
Cost per unit:
- Chicken: $2.50
- Beef: $3.00
- Vegetables: $1.00
Nutritional content per unit:
- Chicken: $250$ calories, $30$g protein, $10$g fat
- Beef: $300$ calories, $20$g protein, $20$g fat
- Vegetables: $50$ calories, $5$g protein, $0$g fat
Minimum daily requirements:
- Calories: $2000$
- Protein: $60$g
- Fat: $50$g

Solution using PuLP:

import pulp

# Define the problem
prob = pulp.LpProblem("Diet_Problem", pulp.LpMinimize)

# Define the decision variables: the amount of each food to consume
food_vars = {
    "Chicken": pulp.LpVariable("Chicken", lowBound=0, cat='Continuous'),
    "Beef": pulp.LpVariable("Beef", lowBound=0, cat='Continuous'),
    "Vegetables": pulp.LpVariable("Vegetables", lowBound=0, cat='Continuous')
}

# Define the cost of each food
costs = {
    "Chicken": 2.50,
    "Beef": 3.00,
    "Vegetables": 1.00
}

# Objective function: Minimize the total cost
prob += pulp.lpSum(food_vars[food] * costs[food] for food in food_vars)

# Define the nutritional content of each food
calories = {
    "Chicken": 250,
    "Beef": 300,
    "Vegetables": 50
}

protein = {
    "Chicken": 30,
    "Beef": 20,
    "Vegetables": 5
}

fat = {
    "Chicken": 10,
    "Beef": 20,
    "Vegetables": 0
}

# Add constraints for minimum daily nutritional requirements
prob += pulp.lpSum(food_vars[food] * calories[food] for food in food_vars) >= 2000, "CaloriesRequirement"
prob += pulp.lpSum(food_vars[food] * protein[food] for food in food_vars) >= 60, "ProteinRequirement"
prob += pulp.lpSum(food_vars[food] * fat[food] for food in food_vars) >= 50, "FatRequirement"

# Solve the problem
prob.solve()

# Print the results
print("Optimal diet:")
for food in food_vars:
    print(f"{food}: {food_vars[food].varValue:.2f} units")

print(f"Total cost: ${pulp.value(prob.objective):.2f}")

Explanation:

Decision Variables: food_vars represent the amount of each food to consume.
Objective Function: Minimize the total cost of the diet.
Constraints: Ensure that the total intake of $calories$, $protein$, and $fat$ meets or exceeds the minimum daily requirements.

Output:

Optimal diet:
Chicken: 0.00 units
Beef: 6.67 units
Vegetables: 0.00 units
Total cost: $20.00

Explanation of the Result:

The optimal solution suggests that the most $cost$-$effective$ way to meet the daily nutritional requirements is to consume $6.67$ units of Beef and no units of Chicken or Vegetables.

The total cost of this diet is $20.00.

Key Points:

Beef-only diet: The solution indicates that Beef alone is sufficient to meet the required minimums for calories, protein, and fat, making it the $cheapest$ option at $3.00 per unit.
No Chicken or Vegetables: Since Beef provides a higher concentration of calories and fat compared to Chicken and Vegetables, the solver determined that only Beef is needed to satisfy the constraints without incurring extra costs.
Total Cost: By consuming $6.67$ units of Beef, the total expenditure on this diet is minimized to $20.00.

This result demonstrates how optimization can find the most economical way to satisfy nutritional needs based on the available options.