What's the best way to handle when a constraint evaluates to infinity in OpenMDAO?

Question

As the question states, I'm wondering the best practice when a constraint or function in general evaluates to infinity.

For example, say we are constraining the factor of safety on a part where the factor is calculated as

safety_factor = np.divide(allowable_force, applied_force)

In the case where the applied force is 0, the safety_factor evaluates to infinity. As the the limit of the safety factor as the applied force approaches 0 is infinity, so mathematically this checks out. Say the safety factor is constrained to be above 2, infinity is greater than 2, therefor the constraint should be satisfied. However in practice this leads to the design variables turning to NaN on the next iteration.

This can be seen in the following simple example code

import openmdao.api as om 
import numpy as np

class A(om.ExplicitComponent):

    def setup(self):
        self.add_input('x')
        self.add_output('y')
        self.add_output('y2')

    def compute(self, inputs, outputs):
        outputs['y'] = inputs['x']**2
        outputs['y2']= np.inf

if __name__ == '__main__':

    prob = om.Problem()
    model = prob.model

    idv = model.add_subsystem('idv', om.IndepVarComp(), promotes=['*'])
    idv.add_output('x')
    model.add_subsystem('A', A(), promotes=['*'])

    model.add_design_var('x')

    model.add_constraint('x', lower=2)
    model.add_constraint('y2', lower=2)
    model.add_objective('y')

    # Setup Optimization
    prob.driver = om.ScipyOptimizeDriver()
    prob.driver.options['optimizer'] = 'SLSQP'
    prob.driver.options['maxiter'] = 5
    prob.driver.options['debug_print'] = ['desvars', 'nl_cons', 'objs', 'totals']
    model.approx_totals()

    prob.setup()
    prob.run_driver()

    print('---')
    print(prob['y'])

Despite y2 always being greater than 2, the design variable x still evaluates to NaN. I've done a quick fix where I just check if the constraint is infinite and replace it with a generic very large float (such as 999999999), however this doesn't seem very pythonic and I'm therfore curious the best practice.

Answer 1

I think the best way to handle it is to avoid it by rearranging the equation. Instead of

con1: (allowable_force / applied_force) > 2

try:

con1: applied_force < 0.5 * allowable_force 

There is no division there, so you don't have any divide by zero.

You definitely want to avoid situations where inf or NaN are expected to regularly appear in your variables because your solvers and optimizers are probably not going to be able to converge or to find the correct search direction.

However, for a component output calculation where it can't be avoided, there have been times where we have had to replace an output with a small value when it was too close to zero. This also introduces a discontinuity in the derivatives, which can also be a problem.