0

検証の手段として OpenMDAO を使用して、 MDO テスト スイートの問題を解決しようとしています。私が直面している問題は、特定の問題のオプティマイザー/アルゴリズムを選択することです。

以下は、 IDF メソッドによって定式化されたスイートからのプロパン燃焼問題のコードです。pyOptSparseDriver で利用可能なSLSQP & NSGA2 を使用して解決しようとしました。しかし、どちらのオプティマイザーも最適ではないソリューションを提供しています。

    from __future__ import print_function, division
    import numpy as np
    import time
    from openmdao.api import *

    class PropaneDis1(Component):
        def __init__(self):
            super(PropaneDis1, self).__init__()

            self.add_param('z1', 1.0)
            self.add_param('y3', np.ones(3))

            self.add_output('y1', val=np.ones(2))

        def solve_nonlinear(self, params, unknowns, resids):

            z1 = params['z1']
            y31 = params['y3'][0]

            y12 = 3.0 - z1
            y11 = z1*y31/y12

            unknowns['y1'] = np.array([y11, y12])

    class PropaneDis2(Component):
        def __init__(self):
            super(PropaneDis2, self).__init__()

            self.add_param('z1', val= 0.0)
            self.add_param('y1', val= np.ones(2))
            self.add_param('y3', val= np.ones(3))

            self.add_output('y2', val=np.ones(2))

        def solve_nonlinear(self, params, unknowns, resids):

            z1 = params['z1']
            y12 = params['y1'][1]
            y33 = params['y3'][2]

            y21 = 0.1*z1*y33/(40.0*y12)
            y22 = 0.1*(z1**2)*y33/(40.0*(y12**2))

            unknowns['y2'] = np.array([y21, y22])

    class PropaneDis3(Component):
        def __init__(self):
            super(PropaneDis3, self).__init__()

            self.add_param('z1', val=0.0)
            self.add_param('y1', val=np.ones(2))
            self.add_param('y2', val=np.ones(2))
            self.add_param('x3', val=np.ones(3))

            self.add_output('y3', val=np.ones(3))

        def solve_nonlinear(self, params, unknowns, resids):

            z1 = params['z1']
            y11 = params['y1'][0]
            x31 = params['x3'][0]
            y12 = params['y1'][1]
            x32 = params['x3'][1]
            x33 = params['x3'][2]
            y21 = params['y2'][0]
            y22 = params['y2'][1]

            y31 = (8.0 - 2*y11 - x32 - x33)/2.0
            y32 = 4*10.0 - 2*x31
            y33 = z1 + y11 + x31 + y12 + x32 + x33 + y21 + y22 + y31 + y32

            unknowns['y3'] = np.array([y31, y32, y33])


    class PropaneIDF(Group):

        def __init__(self):
            super(PropaneIDF, self).__init__()

            self.add('pz1', IndepVarComp('z1', val=1.0), promotes=['*'])
            self.add('px3', IndepVarComp('x3', val=np.array([20.0,1.5,1.0])), promotes=['*'])
            self.add('py1t', IndepVarComp('y1t', val=np.array([1.0,1.0])), promotes=['*'])
            self.add('py2t', IndepVarComp('y2t', val=np.array([1.0,1.0])), promotes=['*'])
            self.add('py3t', IndepVarComp('y3t', val=np.array([1.0,1.0,1.0])), promotes=['*'])

            self.add('d1', PropaneDis1(), promotes=['z1', 'y1'])
            self.add('d2', PropaneDis2(), promotes=['z1', 'y2'])
            self.add('d3', PropaneDis3(), promotes=['z1', 'x3', 'y3'])

            self.connect('y1t', 'd2.y1')
            self.connect('y1t', 'd3.y1')

            self.connect('y2t', 'd3.y2')

            self.connect('y3t', 'd1.y3')
            self.connect('y3t', 'd2.y3')

            self.add('con_cmp1', ExecComp('f2=2*z1 + y1[0] + y1[1] + x3[2] + y2[0] + y3[1] + 2*y2[1] - 10.0'
                     , z1=0.0, y1=np.zeros(2), y2=np.zeros(2), y3=np.zeros(3), x3=np.zeros(3))
                     , promotes=['*'])
            self.add('con_cmp2', ExecComp('f6=abs(y1[0]*y1[1])**0.5 - x3[1]*abs(40.0*z1/y3[2])**0.5'
                     , z1=0.0, y1=np.zeros(2), y3=np.ones(3), x3=np.zeros(3))
                     , promotes=['*'])
            self.add('con_cmp3', ExecComp('f7=abs(z1*y1[0])**0.5 - x3[2]*abs(40.0*y1[1]/y3[2])**0.5'
                     , z1=0.0, y1=np.zeros(2), y3=np.ones(3), x3=np.zeros(3))
                     , promotes=['*'])
            self.add('con_cmp4', ExecComp('f9=z1*(x3[0]**0.5) - y1[1]*y3[1]*(40.0/y3[2])**0.5'
                     , z1=0.0, y1=np.zeros(2), y3=np.zeros(3), x3=np.zeros(3))
                     , promotes=['*'])

            self.add('con_cmp51', ExecComp('con51=y1t[0] - y1[0]', y1t=np.zeros(2), y1=np.zeros(2)), promotes=['*'])
            self.add('con_cmp52', ExecComp('con52=y1t[1] - y1[1]', y1t=np.zeros(2), y1=np.zeros(2)), promotes=['*'])

            self.add('con_cmp61', ExecComp('con61=y2t[0] - y2[0]', y2t=np.zeros(2), y2=np.zeros(2)), promotes=['*'])
            self.add('con_cmp62', ExecComp('con62=y2t[1] - y2[1]', y2t=np.zeros(2), y2=np.zeros(2)), promotes=['*'])

            self.add('con_cmp71', ExecComp('con71=y3t[0] - y3[0]', y3t=np.zeros(3), y3=np.zeros(3)), promotes=['*'])
            self.add('con_cmp72', ExecComp('con72=y3t[1] - y3[1]', y3t=np.zeros(3), y3=np.zeros(3)), promotes=['*'])
            self.add('con_cmp73', ExecComp('con73=y3t[2] - y3[2]', y3t=np.zeros(3), y3=np.zeros(3)), promotes=['*'])

            self.add('obj_cmp', ExecComp('obj=2*z1 + y1[0] + y1[1] + x3[2] + y2[0] + y3[1] + 2*y2[1] - 10.0 +\
                                         abs(y1[0]*y1[1])**0.5 - x3[1]*abs(40.0*z1/y3[2])**0.5+\
                                         abs(z1*y1[0])**0.5 - x3[2]*abs(40.0*y1[1]/y3[2])**0.5+\
                                         z1*(x3[0]**0.5) - y1[1]*y3[1]*(40.0/y3[2])**0.5 '
                     , z1=0.0, y1=np.zeros(2), y2=np.zeros(2), y3=np.zeros(3), x3=np.zeros(3))
                     , promotes=['*'])

            self.fd_options['force_fd']=True
            self.fd_options['step_size']=1.0e-12

    top = Problem()
    top.root = PropaneIDF()

    top.driver = pyOptSparseDriver()
    top.driver.options['optimizer']='SLSQP'

    top.driver.add_desvar('z1', lower=0.0, upper=20.0)
    top.driver.add_desvar('x3', lower=np.array([0.0, 0.0, 0.0]), upper=np.array([20.0, 20.0, 20.0]))
    top.driver.add_desvar('y1t', lower=np.array([0.0, 0.0]), upper=np.array([20.0, 20.0]))
    top.driver.add_desvar('y2t', lower=np.array([0.0, 0.0]), upper=np.array([20.0, 20.0]))
    top.driver.add_desvar('y3t', lower=np.array([0.0, 0.0, 0.0]), upper=np.array([20.0, 20.0, 20.0]))

    top.driver.add_objective('obj')

    top.driver.add_constraint('con51', equals=0.0)
    top.driver.add_constraint('con52', equals=0.0)
    top.driver.add_constraint('con61', equals=0.0)
    top.driver.add_constraint('con62', equals=0.0)
    top.driver.add_constraint('con71', equals=0.0)
    top.driver.add_constraint('con72', equals=0.0)
    top.driver.add_constraint('con73', equals=0.0)

    top.driver.add_constraint('f2', lower=0.0)
    top.driver.add_constraint('f6', lower=0.0)
    top.driver.add_constraint('f7', lower=0.0)
    top.driver.add_constraint('f9', lower=0.0)
    # top.driver.add_constraint('f2', equals=0.0)
    # top.driver.add_constraint('f6', equals=0.0)
    # top.driver.add_constraint('f7', equals=0.0)
    # top.driver.add_constraint('f9', equals=0.0)

    top.setup()
    tt = time.time()

    top.run()

SLSQP のソリューション -

    Optimization using pyOpt_sparse
    ================================================================================

            Objective Function: _objfunc

        Solution:
    --------------------------------------------------------------------------------
        Total Time:                    0.1535
           User Objective Time :       0.0183
           User Sensitivity Time :     0.0806
           Interface Time :            0.0534
           Opt Solver Time:            0.0013
        Calls to Objective Function :      10
        Calls to Sens Function :            7

        Objectives:
            Name        Value        Optimum
            obj        14.3997             0

       Variables (c - continuous, i - integer, d - discrete):
               Name      Type       Value       Lower Bound  Upper Bound
            z1_0      c       2.130827       0.00e+00     2.00e+01
            x3_0      c      20.000000       0.00e+00     2.00e+01
            x3_1      c       0.000000       0.00e+00     2.00e+01
            x3_2      c      -0.000000       0.00e+00     2.00e+01
           y1t_0      c       3.645641       0.00e+00     2.00e+01
           y1t_1      c       0.869179       0.00e+00     2.00e+01
           y2t_0      c       0.180723       0.00e+00     2.00e+01
           y2t_1      c       0.679352       0.00e+00     2.00e+01
           y3t_0      c       1.691004       0.00e+00     2.00e+01
           y3t_1      c      -0.000000       0.00e+00     2.00e+01
           y3t_2      c      20.000000       0.00e+00     2.00e+01

       Constraints (i - inequality, e - equality):
            Name    Type                    Bounds
           con51      i        0.00e+00 <= -0.499951 <= 0.00e+00
           con52      i        0.00e+00 <= 0.000006 <= 0.00e+00
           con61      i        0.00e+00 <= 0.058146 <= 0.00e+00
           con62      i        0.00e+00 <= 0.378850 <= 0.00e+00
           con71      i        0.00e+00 <= 1.336646 <= 0.00e+00
           con72      i        0.00e+00 <= -0.000000 <= 0.00e+00
           con73      i        0.00e+00 <= -7.860081 <= 0.00e+00
             f2       i        0.00e+00 <= -0.000000 <= 1.00e+20
             f6       i        0.00e+00 <= 1.898220 <= 1.00e+20
             f7       i        0.00e+00 <= 2.972127 <= 1.00e+20
             f9       i        0.00e+00 <= 9.529347 <= 1.00e+20

    --------------------------------------------------------------------------------

NSGA2 のソリューション -

    Optimization using pyOpt_sparse
    ================================================================================

            Objective Function: _objfunc

        Solution:
    --------------------------------------------------------------------------------
        Total Time:                  107.5302
           User Objective Time :      62.5869
           User Sensitivity Time :     0.0000
           Interface Time :           41.6987
           Opt Solver Time:            3.2446
        Calls to Objective Function :   99981
        Calls to Sens Function :            0

        Objectives:
            Name        Value        Optimum
            obj        7.67208             0

       Variables (c - continuous, i - integer, d - discrete):
               Name      Type       Value       Lower Bound  Upper Bound
            z1_0      c       2.184051       0.00e+00     2.00e+01
            x3_0      c      20.000000       0.00e+00     2.00e+01
            x3_1      c       0.000000       0.00e+00     2.00e+01
            x3_2      c       0.000000       0.00e+00     2.00e+01
           y1t_0      c       3.757993       0.00e+00     2.00e+01
           y1t_1      c       0.815949       0.00e+00     2.00e+01
           y2t_0      c       0.193249       0.00e+00     2.00e+01
           y2t_1      c       0.746901       0.00e+00     2.00e+01
           y3t_0      c       1.592243       0.00e+00     2.00e+01
           y3t_1      c      -0.000000       0.00e+00     2.00e+01
           y3t_2      c      20.000000       0.00e+00     2.00e+01

       Constraints (i - inequality, e - equality):
            Name    Type                    Bounds
           con51      i        0.00e+00 <= 0.173682 <= 0.00e+00
           con52      i        0.00e+00 <= -0.520903 <= 0.00e+00
           con61      i        0.00e+00 <= -0.122582 <= 0.00e+00
           con62      i        0.00e+00 <= -0.292454 <= 0.00e+00
           con71      i        0.00e+00 <= -0.120713 <= 0.00e+00
           con72      i        0.00e+00 <= 0.000372 <= 0.00e+00
           con73      i        0.00e+00 <= -7.729249 <= 0.00e+00
             f2       i        0.00e+00 <= 0.143730 <= 1.00e+20
             f6       i        0.00e+00 <= 1.641686 <= 1.00e+20
             f7       i        0.00e+00 <= 1.957046 <= 1.00e+20
             f9       i        0.00e+00 <= 3.929613 <= 1.00e+20

    ----------------------------------------------------------------------

制約 f2、f6、f7 & f9 を下限ではなく等式制約 (=0) として問題を実行しようとしたため、問題が定義にないことはわかっています。オプティマを取得しています。

では、他にどのようなオプティマイザーを検討する必要がありますか? または、オプティマイザの問題と混同している他の問題はありますか。

前もって感謝します。

4

1 に答える 1

0

まず、SLSQP の出力を見ると、オプティマイザは実行可能なソリューションをまったく提供できませんでした。具体的には、con51 から con73 までの制約には多数の違反があります。これらは、専門分野コンポーネント間の循環的な依存関係を置き換える互換性の方程式であるため、オプティマイザーは学際的な結合を収束できなかったようです。

私はあなたのコードを取得し、MDF 構成に変換し、非線形ソルバーを Newton に切り替えて、単一の MDA サイクルの収束をどのように処理したかを確認するなど、いくつかのことを試しました (したがって、オプティマイザーも削除しました)。他の設計変数が初期条件のままであると仮定すると、適切な 1e-10 許容誤差に収束するようになりました。

y1 [ 1.  2.]
y2 [ 0.05886036  0.02943018]
y3 [  2.          38.          47.08829054]

y3 の収束値は、IDF で y3t を設計変数として使用するときに指定した上限をはるかに超えていることに注意してください。

top.driver.add_desvar('y3t', lower=np.array([0.0, 0.0, 0.0]),
                             upper=np.array([10.0, 10.0, 10.0]))

プロパンの問題で上限が指定されていないため、どこで上限を取得したかわかりません。さらに重要なことは、あなたのコードに表示されているものと、MDO テスト スイートのリンクで指定されているものを一致させることができませんでした。何らかのエラーがあると思われますが、何も確認できませんでした。

于 2016-05-31T15:21:00.493 に答える