# Example

## Select a system to be optimized and create a designable model using this interface

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

## Instrument the model to be optimization ready

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)

## Create an objective from the metric function and target values

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

## Create an optimizer and give it the objective you defined

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

## Ask the optimizer to find a final design that will optimize the system

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

## Examine the result

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

The L2 norm of the distance between the results.x and exactX is 0.052714038720902294.

You can see that the optimizer did indeed find an optimal solution.

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]

The L2 norm of the distance between the results2.x and exactX is 0.04935654019790089.

You can see that the optimizer did indeed find a better solution, though at the expense of more iterations.