最適化するシステムを選択して、このインターフェースを使用して設計可能モデルを作成します。

この例では、設計可能な4バー機構が使用されています。
################################################################################
# Model definition                                                             #
################################################################################
def fourbar (design):

  """
  The 4 bar is defined with 4 design points A, B, C, D.
 
                                                       coupler (p3)
                                                          /                                                              
  The 4 bar is an RSUR mechanism               B o- - - - - - -o C
    R is at Point A                             /               \
    S is at Point B            ^ Yg            / crank           \follower
    U is at Point C            |              /                   \
    R is at point D            o--->Xg     A o~~~~~~~~~~~~~~~~~~~~~o D

  The two revolute joints at A and D are defined with their Z-axes along the global Z
  The spherical joint at B does not care about orientations
  The universal joint at C is defined as follows:
    1st z-axis, zi,  along global Z 
    2nd z-axis, zj,  perpendicular to zi and the line from B to C
     
  A = (0,0)  B=(40,200)  C=(360,600)  D=(400,0)
  
  The entire model is parameterized in terms of these 4 design points: A, B, C and D.
  Since we are operating in 2D space, this leads to 8 DVs: ax, ay, bx, by, cx, cy, dx & dy.

"""
  # Helper function to create a PART, MARKERS, GRAPHIC, REQUEST
  def make_cylinder (a, b, r, rho, label, oxyz):
    """ Create a cylinder, its CM and graphic given:
        - End points a and b
        - Radius r 
        - density rho
    """
    r2      = r*r
    h2      = (a[0]-b[0])**2 + (a[1]-b[1])**2 + (a[2]-b[2])**2
    h       = math.sqrt(h2)
    mass    = math.pi * r2 * h * rho
    ixx     = iyy = mass * (3 * r2 + h2) / 12
    izz     = mass * r2 / 2
    qploc   = (a + b)/2
    xploc   = qploc + [0,0,1]

    # Create a PART, its CM, a MARKER for CYLINDER reference, a CYLINDER and a REQUEST
    cyl     = Part     (mass=mass, ip=[ixx, iyy, izz, 0, 0, 0], label=label)
    cyl.cm  = Marker   (body=cyl, qp=qploc, zp=b, xp=xploc, label=label + "_CM")  
    cyl.gcm = Marker   (body=cyl, qp=a, zp=b, xp=xploc, label=label + "_Cylinder")
    cyl.gra = Graphics (type="CYLINDER", cm=cyl.gcm, radius=r, length=h, seg=40)
    cyl.req = Request  (type="DISPLACEMENT", i=cyl.cm, comment=label + " CM displacement") 
    return cyl

  # Create a REVOLUTE or SPHERICAL joint
  def make_joint (p1, p2, loc, type):
    xploc = loc + [1,0,0]
    zploc = loc + [0,0,1]
    p1m   = Marker (body=p1, qp=loc, zp=zploc, xp=xploc)
    p2m   = Marker (body=p2, qp=loc, zp=zploc, xp=xploc)
    joint = Joint (type=type,  i=p1m, j=p2m)
    return joint

  model  = Model (output="test")
  units  = Units (mass="KILOGRAM", length="MILLIMETER", time="SECOND", force="NEWTON")
  grav   = Accgrav (jgrav=-9800)
  output = Output (reqsave=True)

  # define the design variables
  ax     = Dv ( label="X coordinate of Point A", b=design["ax"])
  ay     = Dv ( label="Y coordinate of Point A", b=design["ay"])
  bx     = Dv ( label="X coordinate of Point B", b=design["bx"])
  by     = Dv ( label="Y coordinate of Point B", b=design["by"])
  cx     = Dv ( label="X coordinate of Point C", b=design["cx"])
  cy     = Dv ( label="Y coordinate of Point C", b=design["cy"])
  dx     = Dv ( label="X coordinate of Point D", b=design["dx"])
  dy     = Dv ( label="Y coordinate of Point D", b=design["dy"])

  kt, ct, rad, rho = 10, 1, 10, 7.8e-6
  ux, uz           = Vector(1,0,0), Vector(0,0,1) # Global X, # Global Z

  pa     = Point (ax,ay,0)
  pb     = Point (bx,by,0)
  pc     = Point (cx,cy,0)
  pd     = Point (dx,dy,0)

  #Ground
  p0     = Part   (ground=True)
  oxyz   = Marker (body=p0, label="Global CS")

  #Dummy Part  
  p1     = Part  (label="dummy part")
  fj     = Joint (i = Marker(body=p0), j = Marker(body=p1), type="FIXED")
  
  #crank, coupler, follower parts
  p2     = make_cylinder (pa, pb, rad, rho, "Crank")
  p3     = make_cylinder (pb, pc, rad, rho, "Coupler")
  p4     = make_cylinder (pc, pd, rad, rho, "Follower")

  #revolute joint at A, D & spherical at C
  ja     = make_joint (p2, p1, pa, "REVOLUTE")
  jb     = make_joint (p3, p2, pb, "SPHERICAL")
  jd     = make_joint (p1, p4, pd, "REVOLUTE")

  #motion to drive system @A
  j1mot  = Motion   (joint=ja, jtype="ROTATION", dtype="DISPLACEMENT", function="360d*time")
  
  #universal joint at C
  xiloc  = pc + ux                        # Xi along global X
  ziloc  = pc + uz                        # Zi along global Z
  xjloc  = pc + uz                        # Xj along global Z
  zixbc  = uz.cross (pc-pb).normalize ()  # unit vector orthogonal to Zi and BC
  zjloc  = pc+zixbc                       # Zj of universal joint
  p3C    = Marker (body=p3, qp=pc, zp=zjloc, xp=xjloc, label="p3C")
  p4C    = Marker (body=p4, qp=pc, zp=ziloc, xp=xiloc, label="p4C")
  jc     = Joint (type="UNIVERSAL", i=p4C, j=p3C)

  return model

このモデルに対して最適化に向けた準備を行います。

カプラーCM原点のx座標がX-Y平面内の特定の経路をたどるようにします。
  • メトリックを定義します。これは、カプラーCMのx座標を定義する式です。
  • そのターゲット値を定義します。
# Define the “metric” that is to be optimized
expr   = "DX({cm})".format(cm=p3.cm.id)
metric = Variable (label="Coupler DX", function=expr)

# Define the target behavior
splineInput = Variable (label="Spline Independent Variable", function="time")

xval = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 
        1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0]

yval = [201.0000,  78.0588, -35.2336, -102.143, -112.562,  -78.7180,  -23.9943, 36.8618, 
        157.3920, 260.2130, 199.9970,   78.0583, -35.2340, -102.143, -112.562, -78.7178,
        -23.9940,  36.8621, 157.3920,  260.2130, 199.9990]

target = Spline (label="Coupler_DX-Target", x=xval, y=yval)

メトリック関数とターゲット値から目標を作成します。

# Define the cost function
rmsError = RMS_Relative(signal=splineInput, target=target, metric=metric)
rmsError.setup()

最適化エンジンを作成し、定義した目標をこのエンジンに付与します。

# Define the initial design
startDesign = OrderedDict ()
startDesign["ax"]  = -15
startDesign["ay"]  =  15
startDesign["bx"]  =  45 
startDesign["by"]  = 220
startDesign["cx"]  = 345
startDesign["cy"]  = 585
startDesign["dx"]  = 415
startDesign["dy"]  = -15

Create an optimizer
opt = Optimizer (label="Optimize RMS Error", design=startDesign, modelGenerator=fourbar,
                 objective=rmsError, simType="STATICS", simTend=4.0, simDTout=0.01)

システムを最適化する最終設計を見つけるように最適化エンジンに指示します。

result = opt.optimize (maxit=100, accuracy=1.0e-2)

結果を確認します。

# Was the optimization successful?
>>> result.success
True

# How many iterations did the optimizer take?
>>> result.nit
21

# What was the final design?
>>> result.x
[0.01983578586441092, -0.00898945915748043, 39.99727645800783, 199.9829124541905,  
 359.9752404909047, 600.0265028357005, 399.9737539465066, -0.0009001577101342088]
正確な解は次のとおりです:
exactX = [0, 0, 40, 200, 360, 600, 400, 0]

results.xとexactXの間の距離のL2ノルムは0.052714038720902294です。

最適化エンジンが実際に最適な解を見つけたことがわかります。

精度を1e-3に設定して再実行した場合は、26回の反復後に以下が得られます:
>>> result2 = opt.optimize (maxit=100, accuracy=1.0e-3)
>>> result2.x
[-0.00228784659143783, 0.013354616280974424, 40.00015734224888, 200.0138744649529, 359.99195451559245, 599.9837146566314, 400.0200018558606, -0.036468425275888644]

results2.xとexactXの間の距離のL2ノルムは0.04935654019790089です。

反復回数は増えましたが、最適化エンジンは実際に先ほどより優れた解を見つけたことがわかります。