CODEWITHSUNDEEP
optimizationconstraintsopenmdao
Kasia
Kasia

Reputation: 105

Treatment of constraints in SLSQP optimization with openMDAO

With openMDAO, I am using FD derivatives and trying to solve a non-linearly constrained optimization problem with the SLSQP method. Any time the optimizer arrives at a point that violates one of the constraints, it just crashes with the message:

Optimization FAILED. Positive directional derivative for linesearch

For instance, if I intentionally set the initial point to an unfeasible design point, the optimizer performs 1 iteration and exits with the above error (the same happens when I start from a feasible point, but then optimizer arrives at an unfeasible point after a few iterations).

Based on the answer to the question in In OpenMDAO, is there a way to ensure that the constraints are respected before proceeding with a computation?, I'm assuming that raising the AnalysisError exception will not work in my case, is that correct? Is there any other way to prevent the optimizer from going into unfeasible regions or at least backtrack on the linesearch and try a different direction/distance? Or should the SLSQP method be only used when analytic derivatives are available?

Reproducible test case:

import numpy as np
import openmdao.api as om


class d1(om.ExplicitComponent):

        def setup(self):

            # Global design variables
            self.add_input('r', val= [3,3,3])          
            self.add_input('T', val= 20)
            
            # Coupling output
            self.add_output('M', val=0)
            self.add_output('cost', val=0)
            
        def setup_partials(self):
            # Finite difference all partials.
            self.declare_partials('*', '*', method='fd')

        def compute(self, inputs, outputs):
            # define inputs
            r = inputs['r']
            T = inputs['T'][0]
            
            cost =  174.42 * T * (r[0]**2 + 2*r[1]**2 + r[2]**2 + r[0]*r[1] + r[1]*r[2])
            
            M =     456.19 * T * (r[0]**2 + 2*r[1]**2 + r[2]**2 + r[0]*r[1] + r[1]*r[2]) - 599718
                 
            outputs['M'] = M
            outputs['cost'] = cost
            

class MDA(om.Group):

    class ObjCmp(om.ExplicitComponent):

        def setup(self):
            # Global Design Variable
            self.add_input('cost', val=0)

            # Output
            self.add_output('obj', val=0.0)

        def setup_partials(self):
            # Finite difference all partials.
            self.declare_partials('*', '*', method='fd')

        def compute(self, inputs, outputs):
            
            outputs['obj'] = inputs['cost']
           
            
    class ConCmp(om.ExplicitComponent):

        def setup(self):
            # Global Design Variable
            self.add_input('M', val=0)
            
            # Output
            self.add_output('con', val=0.0)

        def setup_partials(self):
            # Finite difference all partials.
            self.declare_partials('*', '*', method='fd')

        def compute(self, inputs, outputs):
            
            # assemble outputs
            outputs['con'] = inputs['M']

    def setup(self):
        
        self.add_subsystem('d1', d1(), promotes_inputs=['r','T'],
                            promotes_outputs=['M','cost'])
                                     
        self.add_subsystem('con_cmp', self.ConCmp(), promotes_inputs=['M'],
                            promotes_outputs=['con'])
                          
        self.add_subsystem('obj_cmp', self.ObjCmp(), promotes_inputs=['cost'],
                           promotes_outputs=['obj'])
                           

# Build the model
prob = om.Problem(model=MDA())
model = prob.model

model.add_design_var('r', lower= [3,3,3], upper= [10,10,10])
model.add_design_var('T', lower= 20, upper= 220)
model.add_objective('obj', scaler=1)
model.add_constraint('con', lower=0)

# Setup the optimization
prob.driver = om.ScipyOptimizeDriver(optimizer='SLSQP', tol=1e-3, disp=True)

prob.setup()
prob.set_solver_print(level=0)
prob.run_driver()

# Printout
print('minimum found at')
print(prob.get_val('T')[0])
print(prob.get_val('r'))

print('constraint')
print(prob.get_val('con')[0])

print('minimum objective')
print(prob.get_val('obj')[0])

Upvotes: 0

Views: 485

Answers (2)

Justin Gray
Justin Gray

Reputation: 5710

Based on your provided test case, the problem here is that your have a really poorly scaled objective and constraint (you also have some very strange coding choices ... which I modified).

Running the OpenMDAO scaling report shows that your objective and constraint values are both around 1e6 in magnitude:

scaling report

This is quite large, and is the source of your problems. A (very rough) rule of thumb is that your objectives and constraints should be around order 1. Thats not hard and fast rule, but is generally a good starting point. Sometimes other scaling will work better, if you have very very larger or small derivatives ... but there are parts of SQP methods that are sensitive to the scaling of objective and constraint values directly. So trying to keep them roughly in the range of 1 is a good idea.

Adding ref=1e6 to both objective and constraints gave enough resolution for the numerical methods to converge the problem:

            Current function value: [0.229372]
            Iterations: 8
            Function evaluations: 8
            Gradient evaluations: 8
Optimization Complete
-----------------------------------
minimum found at
20.00006826587515
[3.61138704 3.         3.61138704]
constraint
197.20821903413162
minimum objective
229371.99547899762

Here is the code I modified (including removing the extra class definitions inside your group that didn't seem to be doing anything):

import numpy as np
import openmdao.api as om


class d1(om.ExplicitComponent):

        def setup(self):

            # Global design variables
            self.add_input('r', val= [3,3,3])          
            self.add_input('T', val= 20)
            
            # Coupling output
            self.add_output('M', val=0)
            self.add_output('cost', val=0)
            
        def setup_partials(self):
            # Finite difference all partials.
            self.declare_partials('*', '*', method='cs')

        def compute(self, inputs, outputs):
            # define inputs
            r = inputs['r']
            T = inputs['T'][0]
            
            cost =  174.42 * T * (r[0]**2 + 2*r[1]**2 + r[2]**2 + r[0]*r[1] + r[1]*r[2])
            
            M =     456.19 * T * (r[0]**2 + 2*r[1]**2 + r[2]**2 + r[0]*r[1] + r[1]*r[2]) - 599718
                 
            outputs['M'] = M
            outputs['cost'] = cost
            

class MDA(om.Group):

    def setup(self):
        
        self.add_subsystem('d1', d1(), promotes_inputs=['r','T'],
                            promotes_outputs=['M','cost'])
                                     
        # self.add_subsystem('con_cmp', self.ConCmp(), promotes_inputs=['M'],
        #                     promotes_outputs=['con'])
                          
        # self.add_subsystem('obj_cmp', self.ObjCmp(), promotes_inputs=['cost'],
        #                    promotes_outputs=['obj'])
                           

# Build the model
prob = om.Problem(model=MDA())
model = prob.model

model.add_design_var('r', lower= [3,3,3], upper= [10,10,10])
model.add_design_var('T', lower= 20, upper= 220)
model.add_objective('cost', ref=1e6)
model.add_constraint('M', lower=0, ref=1e6)

# Setup the optimization
prob.driver = om.ScipyOptimizeDriver(optimizer='SLSQP', tol=1e-3, disp=True)

prob.setup()
prob.set_solver_print(level=0)

prob.set_val('r', 7.65)
prob.run_driver()

# Printout
print('minimum found at')
print(prob.get_val('T')[0])
print(prob.get_val('r'))

print('constraint')
print(prob.get_val('M')[0])

print('minimum objective')
print(prob.get_val('cost')[0])

Upvotes: 1

Justin Gray
Justin Gray

Reputation: 5710

Which SLSQP method are you using? There is one implementation in pyOptSparse and one in ScipyOptimizer. The one in pyoptsparse is older and doesn't respect bounds constraints. The one in Scipy is newer and does. (Yes, its very confusing that they have the same name and share some lineage... but are not the same optimizer any more)

You shouldn't raise an analysis error when you go outside the bounds. If you need strict bounds respecting, I suggest using IPopt from within pyoptsparse (if you can get it to compile) or switching to ScipyOptimizerDriver and its SLSQP implementation.

Based on your question, its not totally clear to me if you're talking about bounds constraints or inequality/equality constraints. If its the latter, then then there isn't any optimizer that would guarantee you remain in a feasible region all the time. Interior point methods like IPopt will stay inside the region much better, but not 100% of the time.

In general, with gradient based optimization its pretty critical that you make your problem smooth and continuous even when its outside the constraint areas. If there are parts of the space that you absolutely can not go into, then you need to make those variables into design variables and use bound constraints. This sometimes requires reformulating your problem formulation a little bit, possibly by adding a kind of compatibility constraint that says "design variable = computed_value". Then you can make sure that the design variable is passed into anything that requires the value to be strictly within a bound, and (hopefully) a converged answer will also satisfy your compatibility constraint.

If you provide some kind of a test case or example, I can amend my answer with a more specific suggestion.

Upvotes: 0

Related Questions