Example

A designable 4-bar mechanism is used for this example.

################################################################################
# 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

We want the x-coordinate of the Coupler CM origin to follow a certain path in the X-Y plane.

• Define the metric, the expression defining the x-coordinate of the Coupler CM.
• Define its target value.

# 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]

The exact answer is:

exactX = [0, 0, 40, 200, 360, 600, 400, 0]

If you set the accuracy to 1e-3 and rerun, after 26 iterations you will obtain:

>>> 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]