# Example

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

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

## Instrument the model to be optimization ready

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

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