Problem

Define a problem using modOpt’s Problem class

This example does not intend to cover all the features of the Problem class. For more details on Problem, please see the API Reference. In this example, we solve a constrained problem given by

\[ \begin{align}\begin{aligned} \underset{x_1, x_2 \in \mathbb{R}}{\text{minimize}} \quad x_1^2 + x_2^2\\\newline \text{subject to} \quad x_1 \geq 0 \newline \quad \quad \quad \quad x_1 + x_2 = 1 \newline \quad \quad \quad \quad x_1 - x_2 \geq 1 \end{aligned}\end{align} \]

We know the solution of this problem is \(x_1=1\), and \(x_2=0\). However, we start from an initial guess of \(x_1=500.0\), and \(x_2=5.0\) for the purposes of this tutorial.

The problem model is written as a subclass of as Problem follows:

import numpy as np
import modopt as mo

# minimize x^2 + y^2 subject to x>=0, x+y=1, x-y>=1.

class QuadraticProblem(mo.Problem):
    def initialize(self):
        self.problem_name = 'quadratic_problem'

    def setup(self):
        self.add_design_variables('x',                              # Design variable name
                                  shape=(2, ),                      # Design variable shape
                                  lower=np.array([0., -np.inf]),    # Design variable lower bounds
                                  vals=np.array([500., 5.]))        # Initial guess
        
        self.add_objective('f')                                     # Objective name

        self.add_constraints('c',                                   # Constraint name
                             shape=(2, ),                           # Constraint shape
                             lower=np.array([1., 1.]),              # Constraint lower bounds
                             upper=np.array([1., np.inf]))          # Constraint upper bounds

    def setup_derivatives(self):
        self.declare_objective_gradient(wrt='x')

        constant_J = np.array([[1.,1.],[1.,-1]])                    # Constant value for the constraint Jacobian
        self.declare_constraint_jacobian(of='c',                    
                                         wrt='x',
                                         vals=constant_J)

    def compute_objective(self, dvs, obj):
        x = dvs['x']
        obj['f'] = np.sum(x**2)

    def compute_objective_gradient(self, dvs, grad):
        grad['x'] = 2 * dvs['x']

    def compute_constraints(self, dvs, cons):
        x   = dvs['x']
        con = cons['c']
        con[0] = x[0] + x[1]
        con[1] = x[0] - x[1]

    # NOTE: compute_constraint_jacobian must be defined even if the Jacobian is constant
    def compute_constraint_jacobian(self, dvs, jac):
        pass

Note that in the code above, we combined all the variables and constraints into single vectors. For examples with multiple variable and constraint vectors, see Examples. This is one of the additional convenience features offered by Problem over ProblemLite, as ProblemLite requires users to combine all variables and constraints into single vectors. Although more verbose than ProblemLite, using the Problem class as shown above can be more beneficial for certain problems, e.g, sparse problems where some of the constraints are only affected by a subset of the variables.

Solve your problem using an optimizer

Once your problem is completely defined within a subclass (here QuadraticProblem) of Problem, create a subclass object to pass to the optimizer. Lastly, import your preferred optimizer from modOpt and solve it, following the standard procedure. Here we will use the SLSQP optimizer from the SciPy library.

# Create a Problem subclass object
prob = QuadraticProblem()

# Setup your preferred optimizer (SLSQP) with the Problem object 
# Pass in the options for your chosen optimizer
optimizer = mo.SLSQP(prob, solver_options={'maxiter':20})

# Check first derivatives at the initial guess, if needed
optimizer.check_first_derivatives(prob.x0)

# Solve your optimization problem
optimizer.solve()

# Print results of optimization
optimizer.print_results()
Setting objective name as "f".

----------------------------------------------------------------------------
Derivative type | Calc norm  | FD norm    | Abs error norm | Rel error norm 
----------------------------------------------------------------------------

Gradient        | 1.0000e+03 | 1.0000e+03 | 1.5473e-05     | 1.5472e-08    
Jacobian        | 2.0000e+00 | 2.0000e+00 | 5.0495e-09     | 2.5248e-09    
----------------------------------------------------------------------------


	Solution from Scipy SLSQP:
	----------------------------------------------------------------------------------------------------
	Problem                  : quadratic_problem
	Solver                   : scipy-slsqp
	Success                  : True
	Message                  : Optimization terminated successfully
	Status                   : 0
	Total time               : 0.006365299224853516
	Objective                : 1.0000000068019972
	Gradient norm            : 2.000000006801997
	Total function evals     : 2
	Total gradient evals     : 2
	Major iterations         : 2
	Total callbacks          : 17
	Reused callbacks         : 0
	obj callbacks            : 5
	grad callbacks           : 3
	hess callbacks           : 0
	con callbacks            : 6
	jac callbacks            : 3
	----------------------------------------------------------------------------------------------------

Scaling API

Please refer to the code snippet below as a guide for scaling the design variables, objective, and constraints independent of their definitions.

Warning

The results provided by the optimizer will always be scaled, while the values from the models will remain unscaled.

def setup(self):
    self.add_design_variables('x',                              # Design variable name
                              shape=(2, ),                      # Design variable shape
                              lower=np.array([0., -np.inf]),    # Design variable lower bounds
                              vals=np.array([500., 5.]),        # Initial guess
                              scaler=2.)                        # Design variable scaler - constant value
                            #   scaler=np.array([1., 2.])       # Design variable scaler - vector value
    
    self.add_objective('f', scaler=5.0)                                     # Objective name and scaler

    self.add_constraints('c',                                   # Constraint name
                         shape=(2, ),                           # Constraint shape
                         lower=np.array([1., 1.]),              # Constraint lower bounds
                         upper=np.array([1., np.inf]),          # Constraint upper bounds
                         scaler=np.array([10., 100.]))          # Constraint scaler - vector value
                        #  scaler=5.0)                          # Constraint scaler - constant value