Optimization failing after very few iterations for nonlinear constraints calculated in a blackbox wrapped in an explicit component

Viewed 39

I have a blackbox solver which is wrapped as explicit component and the objective function and constraints are calculated in the blackbox solver and output. These are taken to a constraint components that has an equality constraint defined such that at any iteration, these constraints are satisifed. I am using finite difference to approximate the partial derivatives. However, I get this SLSQP error "Positive directional derivative for linesearch". From S.O., I understand that this error translates - optimizer could not find a direction to move to and also couldn't verify if the results are minimum. I found that for some iterations derivative is 'None' and it was 'None' at least a few times before it threw this error. Is it because the constraints are calculated in the black box solver? or is it because 'fd' for approximation is not working for non linear constraints? or both? A problem summary is attached for reference. enter image description hereenter image description here

from PowerHarvest import *
from HydroDynamics import *
from SatelliteComms import *
from Propulsion import *
from Constraints import *
from SystemCost import *



class MDA(Group):
    """Multidisciplinary Analysis Group"""

    def __init__(self,  derivative_method='fd', **kwargs):
        super(MDA, self).__init__(**kwargs)
        self.derivative_method = 'fd'
     


    def setup(self):

        cycle = self.add_subsystem('cycle',Group(), promotes = ["*"]) 

        cycle.nonlinear_solver = om.NewtonSolver(solve_subsystems = True)
        cycle.nonlinear_solver.options['atol'] = 1e-6
        cycle.add_subsystem('Hydro', Hydro(),promotes = ["*"]) #This is a blackbox explicit component!
        cycle.add_subsystem('Propulsion_system', Propulsion(),promotes = ["*"])
        cycle.add_subsystem('PowerHarvest_system',PowerHarvest(),promotes = ["*"])
        cycle.add_subsystem("SatelitteComs_system", SatelitteComs(),promotes = ["*"])

        
        cycle.nonlinear_solver.options['atol'] = 1.0e-5
        cycle.nonlinear_solver.options['maxiter'] = 500
        cycle.nonlinear_solver.options['iprint'] = 2

        #Add constraint on the each subsytem if possible 
        #cycle.add_constraint('',om.ex)
       
       
        self.add_subsystem('PowerConstraints_system', PowerConstraints(), promotes=["*"])
        self.add_subsystem('BodyConstraints_system', BodyConstraint(),promotes = ["*"])
       
        self.add_subsystem('SystemCost_system',SystemCost(), promotes = ['*'])

      
        self.add_constraint('A_PV', upper = 100, units = 'm**2')
        #these constraints are output of the blackbox solver!
        self.add_constraint('AreaCon', upper = 0)
        self.add_constraint('massCon',equals = 0)
        self.add_constraint('P_Load', upper = 0) # Solar generates just enough for everything no storing!
        self.add_constraint('DraughtCon', lower = 0.5   )
        self.add_constraint('GMCon', lower = 0.01) #should be positive
        #self.add_constraint('theta', upper = 0.14, lower = 0.1)
        self.add_constraint('Amplitude_Con',upper = -0.1) #amplitude differenc

Added. Run script

import openmdao.api as om
from geom_utils import *
from openmdao.api import Problem, Group, ExplicitComponent,ImplicitComponent, IndepVarComp, ExecComp,\
    NonlinearBlockGS, ScipyOptimizeDriver,NewtonSolver,DirectSolver,ScipyKrylov
import os
import numpy as np
from types import FunctionType
from geom_utils import *
from capytaine.meshes.meshes import Mesh
from pprint import pprint

from PowerHarvest import *
from HydroDynamics import *
from SatelliteComms import *
from Propulsion import *
from Constraints import *
from SystemCost import *
from PEARLMDA import *



if __name__ == '__main__':
    prob = Problem()
    model = prob.model = MDA()

    prob.driver = ScipyOptimizeDriver(optimizer = 'SLSQP')
    # prob.model.nonlinear_solver = om.NonlinearBlockGS()
    #prob.driver.options['optimizer'] = 'COBYLA'
    prob.driver.options['tol'] = 1e-5


    prob.model.add_design_var('Df', lower= 6.0, upper=20.0, units = "m")
    prob.model.add_design_var('tf', lower=1.0, upper=4.0, units = "m")
    #prob.model.add_design_var('submergence', upper = -0.9)
    prob.model.add_design_var('Vs', lower=1, upper=2, units = "m/s") #make sure the lower, upper are according to their units.
    prob.model.add_design_var('ld', lower = 3, upper = 7, units = 'm' )

   
    prob.model.add_objective('cost_per_byte' ) 

   
    newton = om.NewtonSolver(solve_subsystems=True)
    newton.linesearch = om.BoundsEnforceLS()
    prob.model.nonlinear_solver = newton
    prob.model.linear_solver = om.DirectSolver()

    # sqlite file to record the intermediate calculations and derivatives
    r = om.SqliteRecorder("pearl_computations.sql")
    prob.add_recorder(r)
    prob.driver.add_recorder(r)
    prob.driver.recording_options["record_derivatives"] = True
    # Attach recorder to a subsystem
    model.nonlinear_solver.add_recorder(r)

    model.add_recorder(r)
    prob.driver.recording_options["includes"] = ["*"]

# Attach recorder to a solver
    model.nonlinear_solver.add_recorder(r)
    prob.setup()
    prob.set_solver_print(level=2)

# For gradients across the model this will do the finite difference method
    prob.model.approx_totals(method="fd", step=0.1, form="forward", step_calc="abs")   
    prob.run_model()
    prob.run_driver()
    prob.record("final_state")
    print('minimum objective found at')
    print(prob['cost_per_byte'][0])
    print(prob['A_PV'])
    print(f"tf: {prob['tf'][0]}")

    results = dict()
    results['tf'] = prob['tf'][0]
    results['Df'] = prob['Df'][0]
    results['ld'] = prob['ld'][0]
    results['mass'] = prob['Payloadmass'][0]
    results['DraughCon'] = prob['DraughtCon'][0]
    results['AmplitudeCon'] = prob['AmplitudeCon'][0]

    print(results)

Scaling report enter image description here enter image description here

0 Answers
Related