# 例

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

この例では、設計可能な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
- 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です。

>>> 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です。