
[![INFORMS Journal on Computing Logo](https://INFORMSJoC.github.io/logos/INFORMS_Journal_on_Computing_Header.jpg)](https://pubsonline.informs.org/journal/ijoc)

# 2023.0345
# Cone product reformulation for global optimization

This archive is distributed in association with the [INFORMS Journal on
Computing](https://pubsonline.informs.org/journal/ijoc) under the [MIT License](LICENSE).

The software and data in this repository are a snapshot of the software and data
that were used in the research reported on in the paper 
[Cone product reformulation for global optimization](https://doi.org/10.1287/ijoc.2023.0345) by D. Bertsimas, D. de Moor, D. den Hertog, T. Koukouvinos, and J. Zhen. 

## Cite

To cite the contents of this repository, please cite both the paper and this repo, using their respective DOIs.

https://doi.org/10.1287/ijoc.2023.0345

https://doi.org/10.1287/ijoc.2023.0345.cd

Below is the BibTex for citing this snapshot of the repository.

```
@misc{Cone product reformulation for global optimization,
  author =        {D. Bertsimas, D. de Moor, D. den Hertog, T. Koukouvinos and J. Zhen},
  publisher =     {INFORMS Journal on Computing},
  title =         {Cone product reformulation for global optimization},
  year =          {2025},
  doi =           {10.1287/ijoc.2023.0345.cd},
  url =           {https://github.com/INFORMSJoC/2023.0345},
  note =          {Available for download at https://github.com/INFORMSJoC/2023.0345},
}  
```

## Description

The goal of the software is to implement and replicate the computational experiments in the paper [Cone product reformulation for global optimization](https://doi.org/10.1287/ijoc.2023.0345) by D. Bertsimas, D. de Moor, D. den Hertog, T. Koukouvinos, and J. Zhen. 

## Dependencies

For this project, we use

Julia- We use [Julia](https://julialang.org/downloads/oldreleases/) (version 1.5.3) and the Julia package [JuMP.jl](https://github.com/jump-dev/JuMP.jl) (version 0.21.6) to run the scripts.

MOSEK - We use [MOSEK](https://docs.mosek.com/9.2/toolbox.pdf) (version 9.2.45) to solve the convex relaxations by RPT as a (cone) optimization solver.

Gurobi - We use [Gurobi](https://www.gurobi.com/) (version 9.0.2) to solve the generation of the hyperplanes.

BARON - We use [BARON](https://doi.org/10.1007/BF00138693) (version 20.10.16) to solve nonconvex problems as a global solver.

Pyomo - All computations for BARON are implemented using the Python package [pyomo](https://pyomo.readthedocs.io/en/6.4.1/) (version 6.4.1).

## Numerical experiments

This project contains two numerical experiments: Quadratic optimization with exponential cone constraints (Section 6.1 in [Cone product reformulation for global optimization](https://doi.org/10.1287/ijoc.2023.0345)) and Robust palatable diet problem (Section 6.2 in [Cone product reformulation for global optimization](https://doi.org/10.1287/ijoc.2023.0345)). The codes, data, and results for these experiments can be found in the folders `src`, `Data`, and `Results` respectively. 

---

### Quadratic optimization with exponential cone constraints 

#### Algorithm

We have the following functions 

- `calculate_analytic_center_qp_exp(α, β, C, d)`

  **Purpose:**  
  We want the hyperplane generated in the BB step to cut the feasible region intwo two equally spaced sub regions. Hence, we require the hyperplane to pass through the analytic center. The analytic center is defined as a point that maximizes the product of the slacks of the constraints. 

  **Inputs:**
  - `α` – upper bound in `log(t) ≤ α` with `t = ∑ exp(-x_i)`.
  - `β` – upper bound in `exp(x_i) ≤ β` (componentwise).
  - `C` – constraint matrix (size `m × n`); can be empty (`[]`) for no linear constraints.
  - `d` – right-hand side vector (length `m`), aligned with `C`; can be empty (`[]`) for no linear constraints.

  **Output:**
  - `x_ac` – analytic-center point in \(\mathbb{R}^n\).


- `generate_hyperplane_mip_qp_exp(X, α, β, C, d, δ1, use_ac)`

  **Purpose:**  
  Generates a separating hyperplane \(f^\top x = l\) by solving a linear mip with Gurobi that separates the candidate points `X` as much as possible. This is used as the split in the branch-and-bound.

  **Inputs:**
  - `X` – matrix of candidate vectors (size `n × L`), each column a point in \(\mathbb{R}^n\).
  - `α, β, C, d` – same as above (only used if `use_ac == true` to enforce the hyperplane through the analytic center).
  - `δ1` – small positive parameter in the constraint `sum(f) ≥ δ1` to avoid the trivial zero hyperplane.
  - `use_ac` – if `true`, the hyperplane is constrained to satisfy `f' * x_ac == l` where `x_ac` is the analytic center.

  **Outputs:**
  - `res` – flag indicating whether the above function failed and a fallback was used (`true` means failure).
  - `f` – normal vector of the hyperplane.
  - `l` – right-hand-side constant of the hyperplane \(f^\top x = l\).


- `calculate_ub_qp_exp(Y, P, q, r, α, β, C, d)`

  **Purpose:**  
  Computes an upper bound on the objective by evaluating
  \(\,x^\top P x + q^\top x + r\,\) at a collection of candidate points and taking the best value.

  **Inputs:**
  - `Y` – a matrix of candidate points, with each column `Y[:, i]` a vector in \(\mathbb{R}^n\). 
  - `P` – quadratic coefficient matrix (`n × n`).
  - `q` – linear coefficient vector (`n`).
  - `r` – constant term in the objective.
  - `α, β, C, d` – kept for interface consistency; not used in the current implementation.

  **Output:**
  - `best_ub` – smallest objective value over the given candidates.


- `solve_nonconv_quad_exp_new(P, q, r, α, β, C, d, use_lmi, use_new)`

  **Purpose:**  
  Solves a convex conic relaxation of the nonconvex problem using Mosek. The relaxation is written with exponential-cone constraints (and optionally an LMI), providing a **lower bound** and lifted variables that are later used to generate candidate points.

  **Inputs:**
  - `P` – quadratic matrix (`n × n`).
  - `q` – linear vector (`n`).
  - `r` – constant term.
  - `α`, `β` – parameters in the exponential constraints.
  - `C` – optional linear constraint matrix (`m × n`); can be empty (`[]`) if no linear constraints are present.
  - `d` – right-hand side (`m`), aligned with `C`; can be empty (`[]`) if no linear constraints are present.
  - `use_lmi` – if `true`, adds the PSD block constraint on the lifted matrix  
    \([X\; V\; W\; x;\; V'\; Z\; Q\; z;\; W'\; Q'\; T\; t;\; x'\; z'\; t'\; 1] \succeq 0\).
  - `use_new` – switch for formulation variants (kept for compatibility with other code).

  **Outputs:**
  - `x` – relaxed primal vector (candidate solution in \(\mathbb{R}^n\)).
  - `z` – auxiliary vector associated with the `exp(-x)` side.
  - `X` – lifted matrix approximating \(x x^\top\).
  - `V` – lifted matrix associated with cross terms related to `z`.
  - `lb` – lower bound given by the relaxation (`tr(P*X) + q'*x + r`).


- `is_feasible_qp_exp_new(y, α, β)`

  **Purpose:**  
  Checks whether a given vector satisfies the original exponential constraints.

  **Inputs:**
  - `y` – candidate vector in \(\mathbb{R}^n\).
  - `α`, `β` – parameters in the exponential constraints.

  **Constraint check:**
  - \(\log\Big(\sum_i e^{-y_i}\Big) \le α\),
  - \(\sum_i e^{y_i} \le β\).

  **Output:**
  - `is_feas` – `true` if both constraints are satisfied.


- `calculate_cand_vectors_qp_exp_new(x, z, X, V, α, β)`

  **Purpose:**  
  Extracts candidate primal vectors from the relaxation solution \((x, z, X, V)\). The function starts from `x`, and then adds columns `X[:, j]/x[j]` and `V[:, j]/z[j]` that satisfy the exponential constraints.

  **Inputs:**
  - `x` – relaxed primal vector.
  - `z` – auxiliary vector.
  - `X` – lifted matrix associated with `x`.
  - `V` – lifted matrix associated with `z`.
  - `α`, `β` – used for feasibility checks via `is_feasible_qp_exp_new`.

  **Output:**
  - `Y` – list of feasible candidate vectors (each entry is a vector in \(\mathbb{R}^n\)).


- `qp_exp_ub_ipopt(P, q, r, α, β, C, d, x0)`

  **Purpose:**  
  Solves the **original nonconvex problem** with IPOPT starting from an initial point `x0`. This provides a local upper bound on the global optimum (currently commented out in `calculate_ub_qp_exp` but available as an option).

  **Inputs:**
  - `P`, `q`, `r` – objective data.
  - `α`, `β` – exponential constraint parameters.
  - `C`, `d` – optional linear constraints (`Cx ≤ d`), may be empty.
  - `x0` – starting point for IPOPT.

  **Output:**
  - `ub` – local optimum objective value (upper bound).

- `calculate_different_vectors(X, x)`

  **Purpose:**  
  Removes near-duplicate candidate vectors. Starts from a reference vector `x` and concatenates only those candidates that are at least `1e-3` away from `x` in Euclidean norm.

  **Inputs:**
  - `X` – collection of candidate vectors.
  - `x` – reference vector.

  **Output:**
  - `Y` – matrix with first column `x` and subsequent columns distinct from `x`.


- `rpt_branch_and_split_qp_exp_new(P, q, r, α, β, δ, δ1, use_lmi, use_new)`

  **Purpose:**  
  Main **branch-and-bound** algorithm for the QP with exponential constraints.  
  At each node, the algorithm:
  1. Solves the relaxation via `solve_nonconv_quad_exp_new` (lower bound).
  2. Generates candidate points via `calculate_cand_vectors_qp_exp_new`.
  3. Builds a set of distinct candidates via `calculate_different_vectors`.
  4. Computes an upper bound using `calculate_ub_qp_exp`.
  5. Generates a separating hyperplane with `generate_hyperplane_mip_qp_exp` and branches on it.
  6. Repeats until the global gap `UB − LB ≤ δ`, the time limit is reached, or no meaningful split is possible.

  **Inputs:**
  - `P`, `q`, `r` – objective data.
  - `α`, `β` – exponential constraint parameters.
  - `δ` – global optimality tolerance (stop when `UB − LB ≤ δ`).
  - `δ1` – parameter in the hyperplane generation problem (controls nondegeneracy of the normal).
  - `use_lmi` – if `true`, uses the SDP relaxation at every node.
  - `use_newl` – switch for formulation variants (same as in `solve_nonconv_quad_exp_new`).

  **Outputs:**
  - `x_opt` – optimal solutiom.
  - `obj_opt` – objective value at `x_opt` (final upper bound).
  - `gen_hyper` – number of generated hyperplanes (approximate tree depth).
  - `ub_list` – sequence of upper bounds across nodes.
  - `lb_list` – sequence of lower bounds across nodes.
  - `hyp_list_f` – stored hyperplane direction vectors `f`.
  - `hyp_list_l` – stored hyperplane levels `l` (the hyperplane is given by \(f^\top x = l\)).


#### Instances

The parameter values used in all numerical experiments follow exactly the specifications described in the paper. These same values are implemented in the corresponding Julia script located in the src folder.

#### Data

For each instance we provide the objective-function data files `P`, `q`, and `r`, representing `A_0`, `b_0`, and `c` respectively. Each data file contains the exact matrices/vectors/scalar used to generate the experimental results.

---


### Robust palatable diet problem

#### Algorithm

- `rpt_branch_and_split_qp_exp_new(P, q, r, α, β, δ, δ1, use_lmi, use_new)`

  This function implements the RPT-BB algorithm for the **robust nonlinear diet problem** with **quadratic objective**.  

  **Inputs:**
    - `A`, `b`, `h` – palatibility parameters, representing `A`, `b`, and `d` as in the paper.
    - `c` – food-cost vector appearing in the robust constraint.
    - `B` – nutrient matrix (rows = foods, columns = nutrients).
    - `η` – required nutrient intake levels.
    - `ρ` – robustness parameter.
    - `θ` – feasible weight matrix defining the RPT decomposition.
    - `q` – parameter for the power-cone representation of the norm.
    - `λ` – upper bound on total cost used in the robust constraint.
    - `δ` – node termination tolerance (optimality gap).
    - `δ1` – minimum-norm requirement for the hyperplane normal.
    - `use_ac` – whether to initialize the first split with an analytic center.
    - `use_lmi` – whether to include the PSD constraint in the relaxation.
    - `use_new` – toggles between the original and improved RPT decomposition.
    - `rob_nutr` – index of the nutrient used in the robust constraint.

  **Outputs:**
    - `r_opt` – optimal solution.
    - `opt_val` – objective value at `r_opt`.
    - `gen_hyper` – number of generated hyperplanes during branching.
    - `ub_list` – list of upper bounds across explored nodes.
    - `lb_list` – list of lower bounds across explored nodes.
    - `hyp_list_f` – list of hyperplane vectors used for splitting.
    - `hyp_list_l` – list of hyperplane offsets used for splitting.

- Role of Auxiliary Functions

  The remaining functions in this module (such as  
  `solve_rpt_sdp_diet_quad`,  
  `generate_hyperplane_mip_diet`,  
  `calculate_candidate_vectors_diet`,  
  `calculate_diet_ub_quad`,  

  serve **the same purpose as in the previous numerical experiment** (the quadratic problem with exponential‐cone constraints).
  The structure, logic, and interactions between these functions mirror those of the earlier experiment, with the main difference being that the **nonlinear constraints arise from diet-robustness power-cone conditions instead of exponential cones**.

- Meaning of `use_new`

  The flag `use_new` switches between:

    - **RPT formulation without additions** (when `use_new = false`), and  
    - **RPT formulation with additions** (when `use_new = true`).

#### Data

- Palatability parameter estimation.

  The palatability parameters used in the quadratic objective are obtained by solving a regression model. The corresponding code (also provided in `src/Robust_palatable.jl`) is:

  ```julia
  df = CSV.read("/Users/thodoris/Downloads/WFP_dataset.csv", DataFrame; header=true)
  df = select!(df, Not(:Salt))
  X = convert(Matrix, df[:,1:24])
  y = convert(Vector,df[:,25])
  n = size(X)[2]

  m1 = Model(Gurobi.Optimizer)
  @variable(m1, A[1:n,1:n])
  @variable(m1, b[1:n])
  @variable(m1, h)

  @objective(m1, Min, (1/n)*sum(
                      (X[i,:]'*A*X[i,:] + b'*X[i,:] + h - y[i])^2
                      for i in 1:n ) )

  optimize!(m1)

  A_opt, b_opt, h_opt = JuMP.value.(A), JuMP.value.(b), JuMP.value(h)
  ```

  as can be found in the `Robust_palatable.jl` file in the src folder.


- Other model parameters
  All remaining parameters for this numerical experiment are specified in the paper.

- The data for this numerical experiment is given in the folder `Data/Robust palatable diet problem`. It contains three files:
  - `WFP_B_Matrix.csv`: This dataset is from [Maragno et al.](https://pubsonline.informs.org/doi/10.1287/opre.2021.0707), retrieved from the [OptiCL Github repo](https://github.com/hwiberg/OptiCL/tree/main/notebooks/WFP/processed-data). It is used for calculating the objective palatability parameters.
  - `WFP_Costs.csv`: This dataset is from [de Moor et al.](https://www.sciencedirect.com/science/article/pii/S0377221724003333), retrieved from the associated [Zenodo archive](https://zenodo.org/records/8091487). It represents the costs parameters `c_k`.
  - `WFP_dataset.csv`. This dataset is from [Peters et al.](https://pubsonline.informs.org/doi/10.1287/inte.2021.1097), also available on [Zenodo archive](https://zenodo.org/records/8091487). It represents the nutrient coefficient values `beta_kl`.
