search

EmuKit / emukit

A Python-based toolbox of various methods in decision making, uncertainty quantification and statistical emulation: multi-fidelity, experimental design, Bayesian optimisation, Bayesian quadrature, etc.
https://emukit.github.io/emukit/
Apache License 2.0
599 stars 128 forks source link

external objective function evaluation,bayesian optimization #342

Closed AmosJoseph closed 3 years ago

AmosJoseph commented 3 years ago

Hi,I perform bayesian optimization concerning external objective function evaluation. 50 iterations without convergence to a optimum! y=y(xi), i=1,2,3,4,5. [xi : 'domain': (-89,90)] y is known to be betweent about -1 and 0.

The last 10 result of x is: [-21.78456618 -73.00092071 22.18547136 -5.09418206 -66.840349 ] [ 71.06573853 27.49788722 72.22720401 43.18581831 -9.69436757] [-81.03723867 -70.63174164 49.34218906 7.5039747 47.87694802] [ 69.72784349 60.04186887 12.42919718 -5.00636684 84.89919082] [ 74.03608881 28.29894764 46.67773804 41.40603805 -21.24343815] [ 90. -89. 90. -28.94424931 0.90656957] [-74.58716275 -73.2201988 53.57584313 -1.38398485 41.89377634] [ 90. -89. 90. 27.03382007 -42.33288842] [ 90. 90. 90. -49.6513781 3.81251082] [ 90. 90. 20.0884023 -89. 26.63867628]]

The value of y does not converge.The result of y is: [-0.2865 ] [-0.5082 ] [-0.4983 ] [-0.5556 ] [-0.5619 ] [-0.5158 ] [-0.3871 ] [-0.5413 ] [-0.382 ] [-0.2973 ] [-0.4383 ] [-0.43184] [-0.5754 ] [-0.53776] [-0.31811] [-0.50625] [-0.24542] [-0.36656] [-0.2654 ] [-0.1931 ] [-0.23614] [-0.5723 ] [-0.25466] [-0.4048 ] [-0.37171] [-0.36802] [-0.53764] [-0.49476] [-0.49307] [-0.57343] [-0.46464] [-0.31664] [-0.59041] [-0.27351] [-0.52201] [-0.33269] [-0.37579] [-0.24802] [-0.63316] [-0.55878] [-0.55674] [-0.4838 ] [-0.42045] [-0.62791] [-0.56185] [-0.46816] [-0.56181] [-0.6266 ] [-0.25185] [-0.57122] [-0.51804] [-0.2131 ][-0.3723 ] [-0.47639] [-0.53004] [-0.47979] [-0.45373] [-0.62409] [-0.54049] [-0.61303] [-0.62705] [-0.53562]]

Why?Shall I do more iterations?

Many thanks!

Best!

ekalosak commented 3 years ago

Hi, thanks for bringing the question up.

Please post the code you're using to run the optimization - you may be missing a key step, but it's hard to say without seeing the code.

AmosJoseph commented 3 years ago

This is the code,performCAE is the function for exteral evaluation.

-- coding: mbcs --

import sys sys.path.append('C:\Users\Apple1\AppData\Local\Programs\Python\Python37\Lib\site-packages') import numpy as np import GPyOpt from part import from material import from section import from optimization import from assembly import from step import from interaction import from load import from mesh import from job import from sketch import from visualization import from connectorBehavior import *

domain =[{'name': 'var0', 'type': 'continuous', 'domain': (-89,90)}, {'name': 'var1', 'type': 'continuous', 'domain': (-89,90)}, {'name': 'var2', 'type': 'continuous', 'domain': (-89,90)}, {'name': 'var3', 'type': 'continuous', 'domain': (-89,90)}, {'name': 'var4', 'type': 'continuous', 'domain': (-89,90)}]

X_step = np.array([[0,45,-45,90,0], [90,-45,45,0,90], [90,-60,60,30,90], [90,90,0,30,90 ], [90,90,90,90,90], [85,60,30,0,90], [30,60,90,0,30], [90,90,0,90,90], [0,90,90,0,90], [45,-45,30,60,90], [70,-45,0,60,90], [90,45,45,90,90]]) Y_step=np.array([[-0.2865],[-0.5082],[-0.4983],[-0.5556],[-0.5619],[-0.5158], [-0.3871],[-0.5413],[-0.3820],[-0.2973],[-0.4383],[-0.43184]])

inpName = 'Job_linear_bo'

iter_count = 40 current_iter = 0
while current_iter < iter_count: ignored_X = np.array([[0,0,0,0,0]]) bo = GPyOpt.methods.BayesianOptimization(f = None, domain = domain, X = X_step, Y = Y_step, de_duplication = True) x_next = bo.suggest_next_locations(ignored_X = ignored_X)

def performCAE(var):
    var0,var1,var2,var3,var4=var
    mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)
    mdb.models['Model-1'].sketches['__profile__'].CircleByCenterPerimeter(center=(
        0.0, 0.0), point1=(16.25, 13.75))
    del mdb.models['Model-1'].sketches['__profile__']
    mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)
    mdb.models['Model-1'].sketches['__profile__'].CircleByCenterPerimeter(center=(
        0.0, 0.0), point1=(25.0, 10.0))
    mdb.models['Model-1'].sketches['__profile__'].RadialDimension(curve=
        mdb.models['Model-1'].sketches['__profile__'].geometry[2], radius=162.0, 
        textPoint=(49.1602935791016, -11.3555545806885))
    mdb.models['Model-1'].Part(dimensionality=THREE_D, name='Part-1', type=
        DEFORMABLE_BODY)
    mdb.models['Model-1'].parts['Part-1'].BaseShellExtrude(depth=1000.0, sketch=
        mdb.models['Model-1'].sketches['__profile__'])
    del mdb.models['Model-1'].sketches['__profile__']
    mdb.models['Model-1'].parts['Part-1'].PartitionFaceByShortestPath(faces=
        mdb.models['Model-1'].parts['Part-1'].faces.getSequenceFromMask(('[#1 ]', 
        ), ), point1=mdb.models['Model-1'].parts['Part-1'].vertices[0], point2=
        mdb.models['Model-1'].parts['Part-1'].vertices[1])
    mdb.models['Model-1'].Material(name='Material-1')
    mdb.models['Model-1'].materials['Material-1'].Elastic(table=((121000.0, 8600.0, 
        8600.0, 0.253, 0.253, 0.421, 3350.0, 3350.0, 2680.0), ), type=
        ENGINEERING_CONSTANTS)
    mdb.models['Model-1'].parts['Part-1'].Surface(name='Surf-1', side1Faces=
        mdb.models['Model-1'].parts['Part-1'].faces.getSequenceFromMask(('[#1 ]', 
        ), ))
    mdb.models['Model-1'].parts['Part-1'].CompositeLayup(description='', 
        elementType=SHELL, name='CompositeLayup-1', offsetType=MIDDLE_SURFACE, 
        symmetric=True, thicknessAssignment=FROM_SECTION)
    mdb.models['Model-1'].parts['Part-1'].compositeLayups['CompositeLayup-1'].Section(
        integrationRule=SIMPSON, poissonDefinition=DEFAULT, preIntegrate=OFF, 
        temperature=GRADIENT, thicknessType=UNIFORM, useDensity=OFF)
    mdb.models['Model-1'].parts['Part-1'].compositeLayups['CompositeLayup-1'].CompositePly(
        additionalRotationField='', additionalRotationType=ROTATION_NONE, angle=0.0
        , axis=AXIS_3, material='Material-1', numIntPoints=3, orientationType=
        SPECIFY_ORIENT, orientationValue=var0, plyName='Ply-1', region=Region(
        faces=mdb.models['Model-1'].parts['Part-1'].faces.getSequenceFromMask(
        mask=('[#1 ]', ), )), suppressed=False, thickness=0.3, thicknessType=
        SPECIFY_THICKNESS)
    mdb.models['Model-1'].parts['Part-1'].compositeLayups['CompositeLayup-1'].CompositePly(
        additionalRotationField='', additionalRotationType=ROTATION_NONE, angle=0.0
        , axis=AXIS_3, material='Material-1', numIntPoints=3, orientationType=
        SPECIFY_ORIENT, orientationValue=var1, plyName='Ply-2', region=Region(
        faces=mdb.models['Model-1'].parts['Part-1'].faces.getSequenceFromMask(
        mask=('[#1 ]', ), )), suppressed=False, thickness=0.3, thicknessType=
        SPECIFY_THICKNESS)
    mdb.models['Model-1'].parts['Part-1'].compositeLayups['CompositeLayup-1'].CompositePly(
        additionalRotationField='', additionalRotationType=ROTATION_NONE, angle=0.0
        , axis=AXIS_3, material='Material-1', numIntPoints=3, orientationType=
        SPECIFY_ORIENT, orientationValue=var2, plyName='Ply-3', region=Region(
        faces=mdb.models['Model-1'].parts['Part-1'].faces.getSequenceFromMask(
        mask=('[#1 ]', ), )), suppressed=False, thickness=0.3, thicknessType=
        SPECIFY_THICKNESS)
    mdb.models['Model-1'].parts['Part-1'].compositeLayups['CompositeLayup-1'].CompositePly(
        additionalRotationField='', additionalRotationType=ROTATION_NONE, angle=0.0
        , axis=AXIS_3, material='Material-1', numIntPoints=3, orientationType=
        SPECIFY_ORIENT, orientationValue=var3, plyName='Ply-4', region=Region(
        faces=mdb.models['Model-1'].parts['Part-1'].faces.getSequenceFromMask(
        mask=('[#1 ]', ), )), suppressed=False, thickness=0.3, thicknessType=
        SPECIFY_THICKNESS)
    mdb.models['Model-1'].parts['Part-1'].compositeLayups['CompositeLayup-1'].CompositePly(
        additionalRotationField='', additionalRotationType=ROTATION_NONE, angle=0.0
        , axis=AXIS_3, material='Material-1', numIntPoints=3, orientationType=
        SPECIFY_ORIENT, orientationValue=var4, plyName='Ply-5', region=Region(
        faces=mdb.models['Model-1'].parts['Part-1'].faces.getSequenceFromMask(
        mask=('[#1 ]', ), )), suppressed=False, thickness=0.3, thicknessType=
        SPECIFY_THICKNESS)
    mdb.models['Model-1'].parts['Part-1'].compositeLayups['CompositeLayup-1'].ReferenceOrientation(
        additionalRotationField='', additionalRotationType=ROTATION_NONE, angle=0.0
        , axis=AXIS_3, flipNormalDirection=False, flipPrimaryDirection=True, 
        localCsys=None, normalAxisDefinition=SURFACE, normalAxisDirection=AXIS_3, 
        normalAxisRegion=mdb.models['Model-1'].parts['Part-1'].surfaces['Surf-1'], 
        orientationType=DISCRETE, primaryAxisDefinition=EDGE, primaryAxisDirection=
        AXIS_1, primaryAxisRegion=Region(
        edges=mdb.models['Model-1'].parts['Part-1'].edges.getSequenceFromMask(
        mask=('[#2 ]', ), )), stackDirection=STACK_3)
    mdb.models['Model-1'].rootAssembly.DatumCsysByDefault(CARTESIAN)
    mdb.models['Model-1'].rootAssembly.Instance(dependent=ON, name='Part-1-1', 
        part=mdb.models['Model-1'].parts['Part-1'])
    mdb.models['Model-1'].BuckleStep(name='Step-1', numEigen=6, previous='Initial', 
        vectors=12)
    mdb.models['Model-1'].Pressure(createStepName='Step-1', distributionType=
        UNIFORM, field='', magnitude=1.0, name='Load-1', region=Region(
        side1Faces=mdb.models['Model-1'].rootAssembly.instances['Part-1-1'].faces.getSequenceFromMask(
        mask=('[#1 ]', ), )))
    mdb.models['Model-1'].EncastreBC(createStepName='Step-1', localCsys=None, name=
        'BC-1', region=Region(
        edges=mdb.models['Model-1'].rootAssembly.instances['Part-1-1'].edges.getSequenceFromMask(
        mask=('[#4 ]', ), )))
    del mdb.models['Model-1'].boundaryConditions['BC-1']
    mdb.models['Model-1'].EncastreBC(createStepName='Initial', localCsys=None, 
        name='BC-1', region=Region(
        edges=mdb.models['Model-1'].rootAssembly.instances['Part-1-1'].edges.getSequenceFromMask(
        mask=('[#4 ]', ), )))
    mdb.models['Model-1'].DisplacementBC(amplitude=UNSET, createStepName='Initial', 
        distributionType=UNIFORM, fieldName='', localCsys=None, name='BC-2', 
        region=Region(
        edges=mdb.models['Model-1'].rootAssembly.instances['Part-1-1'].edges.getSequenceFromMask(
        mask=('[#1 ]', ), )), u1=SET, u2=SET, u3=UNSET, ur1=SET, ur2=SET, ur3=SET)
    mdb.models['Model-1'].parts['Part-1'].seedPart(deviationFactor=0.1, 
        minSizeFactor=0.1, size=20.0)
    mdb.models['Model-1'].parts['Part-1'].setMeshControls(elemShape=QUAD, regions=
        mdb.models['Model-1'].parts['Part-1'].faces.getSequenceFromMask(('[#1 ]', 
        ), ))
    mdb.models['Model-1'].parts['Part-1'].setElementType(elemTypes=(ElemType(
        elemCode=S4R, elemLibrary=STANDARD, secondOrderAccuracy=OFF, 
        hourglassControl=DEFAULT), ElemType(elemCode=S3, elemLibrary=STANDARD)), 
        regions=(mdb.models['Model-1'].parts['Part-1'].faces.getSequenceFromMask((
        '[#1 ]', ), ), ))
    mdb.models['Model-1'].parts['Part-1'].generateMesh()
    mdb.models['Model-1'].rootAssembly.regenerate()
    mdb.Job(atTime=None, contactPrint=OFF, description='', echoPrint=OFF, 
        explicitPrecision=SINGLE, getMemoryFromAnalysis=True, historyPrint=OFF, 
        memory=90, memoryUnits=PERCENTAGE, model='Model-1', modelPrint=OFF, 
        multiprocessingMode=DEFAULT, name=inpName, nodalOutputPrecision=SINGLE, 
        numCpus=1, numDomains=1, numGPUs=0, queue=None, resultsFormat=ODB, scratch=
        '', type=ANALYSIS, userSubroutine='', waitHours=0, waitMinutes=0)
    #Job and sumbit
    try:
        mdb.jobs[inpName].submit(consistencyChecking=OFF)
        #job.waitForCompletion()
        mdb.jobs[inpName].waitForCompletion()
        odb = session.openOdb(name=inpName+'.odb')
        step = odb.steps.values()[0]
        Descr = step.frames[1].description
        Result = -float(re.split('= ', Descr)[-1])
    except BaseException, e:# in case fail to get result and set the default
        Result = 0.0

    return Result

y_next = performCAE(*x_next.tolist())
X_step = np.vstack((X_step, x_next))
Y_step = np.vstack((Y_step, y_next))

current_iter += 1
inpName = 'Job_linear_bo' + str(current_iter)

print(X_step)
print(Y_step)

AmosJoseph commented 3 years ago

The code is a little long,so it is not very neatly for reading after posted on the webpage.Sorry!

AmosJoseph commented 3 years ago

Hi,main code below,what is wrong?

domain =[{'name': 'var0', 'type': 'continuous', 'domain': (-89,90)}, {'name': 'var1', 'type': 'continuous', 'domain': (-89,90)}, {'name': 'var2', 'type': 'continuous', 'domain': (-89,90)}, {'name': 'var3', 'type': 'continuous', 'domain': (-89,90)}, {'name': 'var4', 'type': 'continuous', 'domain': (-89,90)}]

X_step = np.array([[0,45,-45,90,0], [90,-45,45,0,90], [90,-60,60,30,90], [90,90,0,30,90 ], [90,90,90,90,90], [85,60,30,0,90], [30,60,90,0,30], [90,90,0,90,90], [0,90,90,0,90], [45,-45,30,60,90], [70,-45,0,60,90], [90,45,45,90,90]]) Y_step=np.array([[-0.2865],[-0.5082],[-0.4983],[-0.5556],[-0.5619],[-0.5158], [-0.3871],[-0.5413],[-0.3820],[-0.2973],[-0.4383],[-0.43184]])

iter_count = 40 current_iter = 0 while current_iter < iter_count: ignored_X = np.array([[0,0,0,0,0]]) bo = GPyOpt.methods.BayesianOptimization(f = None, domain = domain, X = X_step, Y = Y_step, de_duplication = True) x_next = bo.suggest_next_locations(ignored_X = ignored_X)

def performCAE(var):
    var0,var1,var2,var3,var4=var
    return Result

y_next = performCAE(*x_next.tolist())
X_step = np.vstack((X_step, x_next))
Y_step = np.vstack((Y_step, y_next))

current_iter += 1
AmosJoseph commented 3 years ago

excuse me,I have a problem about doing Bayesian optimization by Emukit in “Emukit-tutorial-bayesian-optimization-external-objective-evaluation.ipynb”.

target_function, space = forrester_function() space refers to domain of X,but how to define domain of X for external-objective-evaluation without explicit expression?

Best!

AmosJoseph commented 3 years ago

from emukit.core import ContinuousParameter from emukit.core import ParameterSpace from emukit.examples.gp_bayesian_optimization.single_objective_bayesian_optimization import GPBayesianOptimization from emukit.core.loop import UserFunctionResult parameter_space = ParameterSpace([ContinuousParameter('x1', -89, 90), ContinuousParameter('x2', -89, 90), ContinuousParameter('x3', -89, 90), ContinuousParameter('x4', -89, 90), ContinuousParameter('x5', -89, 90)])

X = np.array([[0,45,-45,90,0],
[90,-45,45,0,90],
[90,-60,60,30,90],
[90,90,0,30,90]]) Y = np.array([[0.2865],[0.5082],[0.4983],[0.5556]]) def performCAE(var): var0,var1,var2,var3,var4=var

return Result

inpName = 'Job_linear_emukit_bo'
num_iterations = 40 bo = GPBayesianOptimization(variables_list=[ContinuousParameter('x1', -89, 90), ContinuousParameter('x2', -89, 90), ContinuousParameter('x3', -89, 90), ContinuousParameter('x4', -89, 90), ContinuousParameter('x5', -89, 90)], X=X, Y=Y) results = None

for i in range(num_iterations): X_new = bo.get_next_points(results) Y_new = performCAE(*X_new) results = [UserFunctionResult(X_new[0], Y_new[0])] inpName = 'Job_linear_emukit_bo' + str(i)

X = bo.loop_state.X Y = bo.loop_state.Y print(X) print(Y)

Excuse me,would you like to help me find what is wrong?

apaleyes commented 3 years ago

It is hard to get through these postings! Would be really easier if they were formatted properly. ANyhow, I see you are using external evaluation, calling this function:

def performCAE(var):
    var0,var1,var2,var3,var4=var
    return Result

What's going on here? This isn't complete code is it? Basically if this funciton is all over the place, or very noisy, it would be hard for BO to converge

AmosJoseph commented 3 years ago

performCAE is function for finite element analysis(a solid computational mechanics method). What do you mean by all over the place? It is not noisy to the best of my knowledge.

Below is the complete code : linear_BO_5sym.txt

apaleyes commented 3 years ago

Well, as far as optimizaiton setup goes i don't see any problem. Perhaps the model isn't working well for your data, consider using more advanced modeling, different kernel, etc. Pretty sure the method you are using right now uses RBF kernel, perhaps it isn't a good fit for your model.

AmosJoseph commented 3 years ago

Many thanks for your advice!

Excuse me,how to change modeling method and kernel for bo = GPyOpt.methods.BayesianOptimization(f = None, domain = domain, X = X_step, Y = Y_step, de_duplication = True)？ 6.2 Implementing new acquisitions?Do you have any more example?

apaleyes commented 3 years ago

You need to define a custom model with GPy. In GPyOpt you would then pass it in as model=your_model, in Emukit you would use normal BO interfaces and GPyModelWrapper.

AmosJoseph commented 3 years ago

X=[[0,45,-45,90,0], [90,-45,45,0,90], [90,-60,60,30,90], [90,90,0,30,90 ], [90,90,90,90,90]] Y = [[-0.2865],[-0.5082],[-0.4983],[-0.5556],[-0.5619]] kernel = GPy.kern.Exponential(input_dim=5, variance=1., lengthscale=1.) custom_model = GPy.models.GPRegression(X,Y,kernel) model = custom_model

Pardon,then how to implement external objective function evaluation for this case:

bo = GPyOpt.methods.BayesianOptimization(f = None, domain = domain, X = X_step, Y = Y_step, de_duplication = True)？

apaleyes commented 3 years ago

That's probably not the best place to get support for gpyopt. Here is the tutorial: https://nbviewer.jupyter.org/github/SheffieldML/GPyOpt/blob/devel/manual/GPyOpt_external_objective_evaluation.ipynb

AmosJoseph commented 3 years ago

Actually,I conduct external objective function evaluation according to this tutorial!

But how to use more advanced modeling, different kernel concerning this tutorial?

apaleyes commented 3 years ago

That question is more on GP side, so it's outside of scope of both GPyOpt and Emukit, to be honest. The high level idea is to use kernel that best reflects the knowledge you have about the objective funciton. But I am far from the best person to advise on that. Some initial pointers might be kernel cookbook and GPy docs

As said, this now goes way beyond emukit, so i am closing this issue.

AmosJoseph commented 3 years ago

OK.Many thanks!