Merge pull request #1 from cvxpy/ipopt-interface

Ipopt interface prototype

Author: Transurgeon

cvxpy/reductions/solvers/nlp_solvers/ipopt_nlpif.py

@@ -0,0 +1,84 @@
+"""
+This file is the CVXPY QP extension of the Cardinal Optimizer
+"""
+
+import cvxpy.settings as s
+from cvxpy.reductions.solvers.nlp_solvers.nlp_solver import NLPsolver
+from cvxpy.utilities.citations import CITATION_DICT
+
+
+class IPOPT(NLPsolver):
+    """
+    NLP interface for the IPOPT solver
+    """
+    # Solve capabilities
+    MIP_CAPABLE = True
+
+    # Keyword arguments for the CVXPY interface.
+    INTERFACE_ARGS = ["save_file", "reoptimize"]
+
+    # Map between IPOPT status and CVXPY status
+    STATUS_MAP = {
+                  1: s.OPTIMAL,             # optimal
+                  2: s.INFEASIBLE,          # infeasible
+                  3: s.UNBOUNDED,           # unbounded
+                  4: s.INF_OR_UNB,          # infeasible or unbounded
+                  5: s.SOLVER_ERROR,        # numerical
+                  6: s.USER_LIMIT,          # node limit
+                  7: s.OPTIMAL_INACCURATE,  # imprecise
+                  8: s.USER_LIMIT,          # time out
+                  9: s.SOLVER_ERROR,        # unfinished
+                  10: s.USER_LIMIT          # interrupted
+                 }
+
+    def name(self):
+        """
+        The name of solver.
+        """
+        return 'COPT'
+
+    def import_solver(self):
+        """
+        Imports the solver.
+        """
+        import cyipopt  # noqa F401
+
+    def invert(self, solution, inverse_data):
+        """
+        Returns the solution to the original problem given the inverse_data.
+        """
+        pass
+
+    def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None):
+        """
+        Returns the result of the call to the solver.
+
+        Parameters
+        ----------
+        data : dict
+            Data used by the solver.
+        warm_start : bool
+            Not used.
+        verbose : bool
+            Should the solver print output?
+        solver_opts : dict
+            Additional arguments for the solver.
+        solver_cache: None
+            None
+
+        Returns
+        -------
+        tuple
+            (status, optimal value, primal, equality dual, inequality dual)
+        """
+        pass
+
+    def cite(self, data):
+        """Returns bibtex citation for the solver.
+
+        Parameters
+        ----------
+        data : dict
+            Data generated via an apply call.
+        """
+        return CITATION_DICT["COPT"]

cvxpy/reductions/solvers/nlp_solvers/nlp_solver.py

@@ -0,0 +1,48 @@
+"""
+Copyright 2025, the CVXPY developers
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+"""
+
+from cvxpy.constraints import NonNeg, Zero
+from cvxpy.reductions.solvers.solver import Solver
+
+
+class NLPsolver(Solver):
+    """
+    A non-linear programming (NLP) solver.
+    """
+    # Every QP solver supports Zero and NonNeg constraints.
+    SUPPORTED_CONSTRAINTS = [Zero, NonNeg]
+
+    # Some solvers cannot solve problems that do not have constraints.
+    # For such solvers, REQUIRES_CONSTR should be set to True.
+    REQUIRES_CONSTR = False
+
+    IS_MIP = "IS_MIP"
+
+    def accepts(self, problem):
+        # can accept everything?
+        return True
+
+    def apply(self, problem):
+        """
+        Construct NLP problem data stored in a dictionary.
+        The NLP has the following form
+
+            minimize      f(x)
+            subject to    g^l <= g(x) <= g^u
+                          x^l <= x <= x^u
+        where f and g are non-linear (and possibly non-convex) functions
+        """
+        pass

cvxpy/sandbox/clnlbeam.py

@@ -0,0 +1,54 @@
+
+
+import cvxpy as cp
+
+
+def example_clnlbeam():
+    # Problem parameters
+    N = 1000
+    h = 1 / N
+    alpha = 350
+    
+    # Define variables with bounds
+    t = cp.Variable(N+1)  # -1 <= t <= 1
+    x = cp.Variable(N+1)  # -0.05 <= x <= 0.05
+    u = cp.Variable(N+1)  # unbounded
+    
+    # Define objective function
+    # Minimize: sum of 0.5*h*(u[i+1]^2 + u[i]^2) + 0.5*alpha*h*(cos(t[i+1]) + cos(t[i]))
+    objective_terms = []
+    for i in range(N):  # i from 0 to N-1 (Python 0-indexing)
+        control_term = 0.5 * h * (u[i+1]**2 + u[i]**2)
+        # Note: cos() is non-convex, this may cause solver issues
+        trigonometric_term = 0.5 * alpha * h * (cp.cos(t[i+1]) + cp.cos(t[i]))
+        objective_terms.append(control_term + trigonometric_term)
+    
+    objective = cp.Minimize(cp.sum(objective_terms))
+    
+    # Define constraints
+    constraints = []
+    
+    # Variable bounds
+    constraints.extend([
+        t >= -1,
+        t <= 1,
+        x >= -0.05,
+        x <= 0.05
+    ])
+    
+    # Dynamics constraints
+    for i in range(N):  # i from 0 to N-1
+        # x[i+1] - x[i] - 0.5*h*(sin(t[i+1]) + sin(t[i])) == 0
+        # Note: sin() is also non-convex
+        position_constraint = (x[i+1] - x[i] - 
+                             0.5 * h * (cp.sin(t[i+1]) + cp.sin(t[i])) == 0)
+        constraints.append(position_constraint)
+        
+        # t[i+1] - t[i] - 0.5*h*u[i+1] - 0.5*h*u[i] == 0
+        angle_constraint = (t[i+1] - t[i] - 
+                          0.5 * h * u[i+1] - 0.5 * h * u[i] == 0)
+        constraints.append(angle_constraint)
+    
+    # Create and solve the problem
+    problem = cp.Problem(objective, constraints)
+    return problem

cvxpy/sandbox/mle_estimation.py

@@ -0,0 +1,55 @@
+import cyipopt
+import numpy as np
+from reduction_classes import HS071, Bounds_Getter
+
+import cvxpy as cp
+
+
+def example_mle():
+    # Generate data (matching Julia's Random.seed!(1234))
+    n = 1000
+    np.random.seed(1234)
+    data = np.random.randn(n)
+    print(data.var(), data.mean())
+    
+    # Define variables - matching Julia exactly
+    mu = cp.Variable(1, name="mu")  # mean parameter
+    sigma = cp.Variable(1, name="sigma")  # standard deviation, σ >= 0
+    
+    constraints = [mu == sigma**2]
+    
+    # Calculate the objective exactly as in Julia
+    # n / 2 * log(1 / (2 * π * σ^2)) - sum((data[i] - μ)^2 for i in 1:n) / (2 * σ^2)
+    # = n/2 * log(1/(2*π)) - n*log(σ) - sum((data[i] - μ)^2) / (2*σ^2)
+    
+    # Sum of squared residuals
+    residual_sum = cp.sum_squares(data - mu)
+    
+    # The complete objective (including the constant term that Julia includes)
+    log_likelihood = (n / 2) * cp.log(1 / (2 * np.pi * (sigma)**2)) - residual_sum/(2 * (sigma)**2)
+    
+    objective = cp.Minimize(-log_likelihood)
+    
+    # Create problem with constraints
+    problem = cp.Problem(objective, constraints)
+    return problem
+
+bounds = Bounds_Getter(example_mle())
+x0 = [1.0, 0.0]
+
+nlp = cyipopt.Problem(
+   n=len(x0),
+   m=len(bounds.cl),
+   problem_obj=HS071(bounds.new_problem),
+   lb=[1e-6, None],
+   ub=None,
+   cl=bounds.cl,
+   cu=bounds.cu,
+)
+
+nlp.add_option('mu_strategy', 'adaptive')
+nlp.add_option('tol', 1e-7)
+nlp.add_option('hessian_approximation', "limited-memory")
+
+x, info = nlp.solve(x0)
+print(x)

cvxpy/sandbox/portfolio_opt.py

@@ -0,0 +1,54 @@
+import cyipopt
+import numpy as np
+import pandas as pd
+from reduction_classes import HS071, Bounds_Getter
+
+import cvxpy as cp
+
+
+def portfolio_example():
+        # Data
+        df = pd.DataFrame({
+        'IBM': [93.043, 84.585, 111.453, 99.525, 95.819, 114.708, 111.515,
+                113.211, 104.942, 99.827, 91.607, 107.937, 115.590],
+        'WMT': [51.826, 52.823, 56.477, 49.805, 50.287, 51.521, 51.531,
+                48.664, 55.744, 47.916, 49.438, 51.336, 55.081],
+        'SEHI': [1.063, 0.938, 1.000, 0.938, 1.438, 1.700, 2.540, 2.390,
+                3.120, 2.980, 1.900, 1.750, 1.800]
+        })
+
+        # Compute returns
+        returns = df.pct_change().dropna().values
+        r = np.mean(returns, axis=0)
+        Q = np.cov(returns.T)
+
+        # Single-objective optimization
+        x = cp.Variable(3)
+        variance = cp.quad_form(x, Q)
+        expected_return = r @ x
+
+        prob = cp.Problem(
+        cp.Minimize(variance),
+        [cp.sum(x) <= 1000, expected_return >= 50]
+        )
+        return prob
+
+bounds = Bounds_Getter(portfolio_example())
+x0 = [10.0, 10.0, 10.0]
+
+nlp = cyipopt.Problem(
+   n=len(x0),
+   m=len(bounds.cl),
+   problem_obj=HS071(bounds.new_problem),
+   lb=[0.0, 0.0, 0.0],
+   ub=None,
+   cl=bounds.cl,
+   cu=bounds.cu,
+)
+
+nlp.add_option('mu_strategy', 'adaptive')
+nlp.add_option('tol', 1e-7)
+nlp.add_option('hessian_approximation', "limited-memory")
+
+x, info = nlp.solve(x0)
+print(x)

cvxpy/sandbox/qcp_example.py

@@ -0,0 +1,48 @@
+
+import cyipopt
+import numpy as np
+from reduction_classes import HS071, Bounds_Getter
+
+import cvxpy as cp
+
+
+def example_qcp():
+    # Define variables
+    x = cp.Variable(1)
+    y = cp.Variable(1, bounds=[0, np.inf])  # y >= 0
+    z = cp.Variable(1, bounds=[0, np.inf])  # z >= 0
+    
+    # Define objective (maximize x)
+    objective = cp.Minimize(-x)
+    
+    # Define constraints - exact same as JuMP model
+    constraints = [
+        x + y + z == 1,                # Linear equality constraint
+        x**2 + y**2 - z**2 <= 0,      # Quadratic constraint: x*x + y*y - z*z <= 0
+        x**2 - y*z <= 0               # Quadratic constraint: x*x - y*z <= 0
+    ]
+    
+    # Create and solve problem
+    problem = cp.Problem(objective, constraints)
+    return problem
+
+bounds = Bounds_Getter(example_qcp())
+x0 = [0.2, 0.2, 0.2]
+
+nlp = cyipopt.Problem(
+   n=len(x0),
+   m=len(bounds.cl),
+   problem_obj=HS071(bounds.new_problem),
+   lb=bounds.lb,
+   ub=bounds.ub,
+   cl=bounds.cl,
+   cu=bounds.cu,
+)
+
+nlp.add_option('mu_strategy', 'adaptive')
+nlp.add_option('tol', 1e-7)
+nlp.add_option('hessian_approximation', "limited-memory")
+nlp.add_option('print_level', 7)  # Increase for more detailed output
+
+x, info = nlp.solve(x0)
+print(x)