MDO Test suit: Golinski’s Speed Reducer problem using OpenMDAO

Question

Recently, I have started working in the field of multidisciplinary design optimization. I am using OpenMDAO framework for weight optimization of Golinski's speed reducer in the MDO test suit. I want to apply MDF architecture for this problem. I am referring the paper "Benchmarking Multidisciplinary Design Optimization Algorithms" by Tedford and Martins for problem formulation and decomposition. Where they have decomposed this problem into three disciplines and their individual constraints.

While coding, I referred Seller problem from OpenMDAO documentation. I have made three disciplines and a group (speed_mda()) to implement multidisciplinary analysis. I have added objective function and constraints in the speed_mda() group as a subsystem. I have made their connections with discipline outputs (coupled variables). But I have not applied these constraints on individual discipline (Actually, I don't know how to do it). So I have applied all of them on top level group. I am getting output 2713.678 by violating some of the constraints.

Here is my code:

# Discipline 1
class speed_1(om.ExplicitComponent):
    def setup(self):
        self.add_input('z1', val = 0)
        self.add_input('z2', val = 0)
        self.add_output('y1', val = 1)
    def setup_partials(self):
        # Finite difference all partials.
        self.declare_partials('*', '*', method='fd')
    def compute(self,inputs,outputs):
        outputs['y1'] = max(27/(inputs['z1']**2*(inputs['z2'])),397.5/(inputs['z1']**2*inputs['z2']**2), 5*inputs['z1'], 2.6)
# Discipline 2
class speed_2(om.ExplicitComponent):
    def setup(self):
        self.add_input('z1', val = 0)
        self.add_input('z2', val = 0)
        self.add_input('x21', val = 0)
        self.add_output('y2', val = 0)
    def setup_partials(self):
        self.declare_partials('*', '*', method='fd')
    def compute(self,inputs,outputs):
        outputs['y2'] = max((1.93*inputs['x21']**3/(inputs['z1']*inputs['z2']))**0.25,     1/(0.5*(((1.69*10**7)**2)*inputs['x21']**2/(inputs['z1']**2*inputs['z2']**2) + 745)**0.5)**(0.3333), 2.9)
# Discipline 3
class speed_3(om.ExplicitComponent):
    def setup(self):
        self.add_input('z1', val = 0)
        self.add_input('z2', val = 0)
        self.add_input('x31', val = 0)
        self.add_output('y3', val = 0)
    def setup_partials(self):
        self.declare_partials('*', '*', method='fd')
    def compute(self,inputs,outputs):
        outputs['y3'] = max((1.93*inputs['x31']**3/(inputs['z1']*inputs['z2']))**0.25, 1/(85*(((1.69*10**7)**2)*inputs['x31']**2/(inputs['z1']**2*inputs['z2']**2) + 1.575*(10**8))**0.5)**(0.3333), 5)

class speed_mda(om.Group):
    def setup(self):
        # Adding all discipline to MDA
        cycle = self.add_subsystem('cycle',om.Group(),promotes_inputs=['z1', 'z2', 'x21', 'x31'])
        cycle.add_subsystem('d1', speed_1(), promotes_inputs = ['z1', 'z2'])
        cycle.add_subsystem('d2', speed_2(), promotes_inputs=['z1','z2', 'x21'])
        cycle.add_subsystem('d3', speed_3(), promotes_inputs=['z1','z2', 'x31'])
    
        # No need of connections for the discipline
        cycle.set_input_defaults('x21', 7.8)
        cycle.set_input_defaults('x31', 8.3)
        cycle.set_input_defaults('z1', 0.75)
        cycle.set_input_defaults('z2', 22.0)
    
        # Add solver to MDA: Nonlinear Block Gauss Seidel is a gradient free solver
        cycle.nonlinear_solver = om.NonlinearBlockGS()
    
        # Adding obj. function and constraints as a subsystem
        self.add_subsystem('obj_fun',om.ExecComp('obj = (0.7854*y1*z1**2)*(3.3333*z2**2+14.933*z2-43.0934) - 1.5079*y1*(y2**2+y3**2)+7.477*(y2**3+y3**3)+0.7854*(x21*y2**2+x31*y3**2)', z1=0.0,z2=0.0,x21=0.0,x31=0.0), promotes=['x21','x31','z2','z1','obj'])
        self.add_subsystem('con1',om.ExecComp('c1 = z1*z2 - 40.0'), promotes=['c1']) # Global
        self.add_subsystem('con10',om.ExecComp('c10 = y1 - 12.0*z1'), promotes=['c10']) # 1
        self.add_subsystem('con11',om.ExecComp('c11 = y1 - 3.6'), promotes=['c11'])
        self.add_subsystem('con12',om.ExecComp('c12 = y2 - 3.9'), promotes=['c12'])  # 2
        self.add_subsystem('con13',om.ExecComp('c13 = 2.85*y2 - x21'), promotes=['c13'])
        self.add_subsystem('con14',om.ExecComp('c14 = y3 - 5.5'), promotes=['c14'])  # 3
        self.add_subsystem('con15',om.ExecComp('c15 = 2.09*y3 - x31'), promotes=['c15'])
    
        # Connect outputs from MDA (coupled variables) to obj. function and constraints
        self.connect('cycle.d1.y1',['obj_fun.y1','con10.y1','con11.y1'])
        self.connect('cycle.d2.y2',['obj_fun.y2','con12.y2','con13.y2'])
        self.connect('cycle.d3.y3',['obj_fun.y3','con14.y3','con15.y3'])
    
    
# Form topmost group (Problem) and add above MDF model to it
prob = om.Problem()
model = prob.model = speed_mda()
prob.driver = om.ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'
prob.driver.options['tol'] = 1e-9
prob.driver.options['disp'] = True
model.add_design_var('x21', lower=7.3, upper=8.3)
model.add_design_var('x31', lower=7.3, upper=8.3)
model.add_design_var('z1', lower=0.7, upper=0.8)
model.add_design_var('z2', lower=17.0, upper=28.0)
model.add_objective('obj')
model.add_constraint('c1', upper = 0)
model.add_constraint('c10',upper = 0)
model.add_constraint('c11',upper = 0)
model.add_constraint('c12',upper = 0)
model.add_constraint('c13',upper = 0)
model.add_constraint('c14',upper = 0)
model.add_constraint('c15',upper = 0)

prob.model.approx_totals()

prob.setup()
prob.set_solver_print(level=0)
prob.set_val('x21', 7.8)
prob.set_val('x31', 8.3)
prob.set_val('z1', 0.75)
prob.set_val('z2', 22.0)

prob.run_model()
prob.run_driver()

print('minimum found at')
print((prob.get_val('z1')[0],prob.get_val('z2')[0],prob.get_val('x21')[0],prob.get_val('x31')[0],prob.get_val('cycle.d1.y1')[0], prob.get_val('cycle.d2.y2')[0],prob.get_val('cycle.d3.y3')[0]))
print('minumum objective')
print(prob.get_val('obj')[0])

I am getting following output:

Positive directional derivative for linesearch    (Exit mode 8)
            Current function value: [2713.67806668]
            Iterations: 8
            Function evaluations: 4
            Gradient evaluations: 4
Optimization FAILED.
Positive directional derivative for linesearch
-----------------------------------
minimum found at
(0.7, 17.00007920093152, 7.300007915120889, 7.300015825246984, 3.5, 2.9, 5.0)
minumum objective
2713.6780666813797

Here, objective value is far less than the actual and it says optimization has failed. I looked for the above error but my initial guess is inside the bound. My all problem constraints are not satisfied. I also tried to apply constraints on individual disciplines but I couldn't do it. I don't know what is the actual problem, probably I am making some basic conceptual mistake. Can please anyone help me with this.

Answer 1

There is something strange about the way the speed reducer problem is posed in that paper. The use of max function is not technically differentiable at all, which makes this an less-than-ideal way to implement a problem for use in gradient based optimization

Also, that formulation looks really different from other descriptions of the speed-reducer problem I've seen in other papers. This formulation traces back to original work on an MDO test problem suite, where problems that were not really multidisciplinary were broken up into separate blocks and additional constraints were added to ensure compatibility. In the case of this paper, I think the changes to the problem formulation resulted in some less than ideal problem structure. I recommend you look toward a more well posed formulation such as this one

Regardless, when I set the given inputs from the states optimum of their paper I don't get the optimum value that they reported. I get a lower number, so there must be something subtly different about your code. Somehow your code is returning lower values, so check your equations carefully.

I ran your code on OpenMDAO V3.8 and got the following:

/Users/jsgray/work/packages/OpenMDAO/openmdao/core/total_jac.py:1713: UserWarning:Constraints or objectives ['con1.c1', 'con12.c12', 'con13.c13', 'con14.c14', 'con15.c15'] cannot be impacted by the design variables of the problem.
Positive directional derivative for linesearch    (Exit mode 8)
            Current function value: [2713.66402046]
            Iterations: 9
            Function evaluations: 5
            Gradient evaluations: 5
Optimization FAILED.
Positive directional derivative for linesearch
-----------------------------------
minimum found at
(0.7, 17.000000001268262, 7.3, 7.3, 3.5, 2.9, 5.0)
minumum objective
2713.6640204584155

So I see the same value you do, but I also got a helpful additional warning (newly added in V3.8) about constraints that are not affected by any of the design variables. When I commented out those constraints, the result changed to

Optimization terminated successfully    (Exit mode 0)
            Current function value: [2713.66402024]
            Iterations: 10
            Function evaluations: 6
            Gradient evaluations: 6
Optimization Complete
-----------------------------------
minimum found at
(0.7, 17.0, 7.3, 7.3, 3.5, 2.9, 5.0)
minumum objective
2713.6640202393

Which is the same answer as before, but sans the scary warnings from the optimizer. So the error you were seeing was due to having a large number of constraints that, while inherently satisfied, were not actually controllable by the optimizer. This causes rows of all zeros to show up in the Jacobian and hence makes the optimization problem singular. While SLSQP was able to work around the singularity, it caused enough numerical headaches to cause it to throw the warning.