Force: Penalty

Model Element

Description

Force_Penalty defines a generalized force whose purpose is to maintain a “soft” constraint in the system. This means that when the constraint is satisfied, there is no penalty force. However, if the constraint is violated, a penalty force that tries to reduce the violation is generated. The generalized force acts on all the coordinates involved in the definition of the constraint.

Format

<Force_Penalty

id = "integer"

[label = "string"]

penalty = "double"

[penalty1 = "double"]

[[

unilateral = "False" | "True"

smoothing_factor = "real"

]]

{

type = "EXPRESSION"

expr = "motionsolve_expression"

type = "USERSUB"

usrsub_dll_name = "valid_path_name"

usrsub_fnc_name = "custom_fnc_name"

usrsub_param_string = "USER(par_1, ..., par_n)"

type = "USERSUB"

script_name = "valid_path_name"

interpreter = "PYTHON" | "MATLAB"

usrsub_fnc_name = "custom_fnc_name"

usrsub_param_string = "USER(par_1, ..., par_n)"

}

Attributes

Specifies the element identification number (integer>0). This number is unique among all Reference_Variable elements.

label

This attribute describes the name of the Reference_Variable element. The description is primarily used to make the input more readable.

penalty

Specifies a penalty factor to be used in calculating the restoring force that is used to enforce that the algebraic constraint is always zero. Penalty > 0. See Comment 2 for more details on how this attribute is used.

penalty1

Specifies a second penalty factor to be used in calculating the restoring force that is used to enforce that the time derivative of the algebraic constraint. The default value for Penalty1 is zero. When specified, Penalty1 ≥ 0. See Comment 3 for more details on how this attribute is used.

unilateral

A Boolean attribute. Select one from TRUE and FALSE. Indicates whether an equality constraint, or an inequality constraint is being specified. Unilateral is an optional attribute. If defaults to FALSE when not specified. For more information about equality and inequality constraints please see Comments 1 and 4.

smoothing_factor

The penalty force due to a unilateral constraint is ramped up in a smooth manner using a STEP function. This specifies the x-value at which the smoothing is completed. Smoothing_Factor > 0. For more information about smoothing, please see Comment 4.

type

Select from EXPRESSION and USERSUB. Specifies how the constraint governing the penalty force is defined.

•	EXPRESSION specifies that the constraint is defined in a MotionSolve expression that can be evaluated at run-time.

•	USERSUB indicates that the algebraic constraint is specified in a user-defined subroutine.

expr

Specifies the MotionSolve expression that defines the algebraic constraint. Use this attribute only when type= EXPRESSION. Any valid run-time MotionSolve expression can be provided as input.

usrsub_dll_name

Specifies the path and name of the DLL or shared library containing the user subroutine. MotionSolve uses this information to load the user subroutine in the DLL at run time. Use this keyword only when type = USERSUB.

usrsub_fnc_name

This parameter allows you to specify the name for the user subroutine. The default name, PFOSUB, is used when the attribute is not specified. Use this keyword only when type = USERSUB.

usrsub_param_string

The list of parameters that are passed from the data file to the user-defined PFOSUB. Use this keyword only when type = USERSUB.

script_name

Specifies the path and name of the user written script that contains the function specified by usrsub_fnc_name. Use this keyword only when type = USERSUB.

interpreter

Specifies the interpreted language that the user script is written in. Valid choices are MATLAB or PYTHON. Use this keyword only when type = USERSUB.

Comments

1.	Force_Penalty may be used to define unilateral or inequality constraints. Let q be any set of values dependent on system displacements.

•	Equality or bilateral constraints have the form g(q,t) = 0

•	Inequality or unilateral constraints have the form g(q,t) ≤ 0

•	In both cases, a penalty force is generated only when the constraint is violated.

o	For equality constraints when g(q,t) ≠ 0

o	For unilateral constraints when g(q,t) > 0

2.	Assume that a constraint g(q,t)=0 is to be maintained in the model.

One approach to solve the problem is to define a user-defined constraint and to provide the constraint g(q,t)=0 to MotionSolve. The constraint is enforced by an unknown Lagrange multiplier λ that is solved with the rest of the system variables.

•	λ may be a force or a torque.

•	λ acts along the gradient of the constraint, .

•	The generalized force acting on the system is .

For more information on the above-mentioned approach, see Constraint_General.

Force_Penalty uses a slightly different method to accomplish the same. The constraint is “soft”, in other words, it is not strictly enforced. Violations of the constraint are allowed, but a restoring force tending to decrease the constraint violation is internally generated.

•	The penalty force may be a force or a torque.

•	The penalty force has the magnitude, F = -penalty*g(q,t).

•	The penalty force acts along the gradient of the constraint, .

•	The generalized force acting on the system is .

•	Generally speaking, a large penalty ensures a small constraint violation.

There is a close correspondence between a hard constraint with a Lagrange Multiplier and a soft constraint with a penalty force.

3.	If g(q,t)=0, then it automatically implies ġ(q̇,,q,t)=0, where ġ(...) = . q̇ are the velocities corresponding to q.

To enforce ġ(q̇,,q,t)=0, a second component –penalty1*ġ(q̇,,q,t) must be added to the penalty force computed earlier. To enforce both g(q,t)=0 and ġ(q̇,,q,t)=0, the following strategy is used:

•	The penalty force has the magnitude, F = -penaltyg(q,t) – penalty1ġ(q̇,,q,t).

•	The penalty force acts along the gradient of the constraint, .

•	The generalized force acting on the system is .

When penalty1 is specified, MotionSolve uses the above equations to compute the penalty force and the generalized forces on the equations of motion.

4.	A smoothing function is used to smooth the equations in 3 as follows:

•	p = penalty * STEP(g(q),0,0,smoothing_factor,1)

•	p1 = penalty * STEP(g(q),0,0,smoothing_factor,1)

•	F = - pg(q,t) – p1ġ(q̇,,q,t) )

•	The STEP function is used to gradually ramp up the penalty factors as the violations increase. This is shown in Figure 1 below. The transition will become steeper as smoothing factor becomes smaller.

force_penalty_fig1

Figure 1: The smoothing factor for unilateral constraints

5.	Three approaches are available for defining the function g(y,t).

•	As a function expression in the XML file: This is probably the simplest approach, and g(y,t) is defined as an expression composed from the built-in functions and operators in the function expression language

•	As compiled user-written subroutine, provided in a DLL. Any language may be used. MotionSolve provides an interface to Fortran, C and C++. The default name for the user-written subroutine is PFOSUB. However, you are allowed to give the function any name and tell MotionSolve what it is.

•

As a script written in the Python or MATLAB language: Though this approach is slightly less efficient than writing the user-written subroutine in a compiled language, it provides two significant benefits: (A) It provides for a vastly more flexible and powerful environment to compose the the function g(y,t) than function expressions, and, (B) There is no need for compilers and linkers as is the case when Fortran, C or C++ are used.

6.	Excessively large values for penalty and penalty1 can cause numerical difficulties. Think of penalty as being similar to a stiffness and penalty1 to damping.

7.	Very small values of smoothing_factor may also cause numerical difficulties. A good rule of thumb is to make the smoothing factor a fraction of the maximum penetration expected.

8.	Constraints of the form h(q,t) ≥ 0 may be implemented as follows:

•	Define g(q,t) = -h(q,t)

•	Impose the penalty constraint on g(q,t) ≤ 0

Examples

Example-1:

Assume a bead, modeled as a particle, is sliding on a catenary fixed in the global x-y plane as shown in Figure 2 below. The path of the bead is shown as a red circle, and the catenary as a blue line.

force_penalty_fig2

Figure 2: A bead sliding on a catenary

The constraint is represented as: g(x,y) = y – cosh(x) = 0.

The time derivative of the constraint is: ġ(ẋ,ẏ,x,y) = ẏ - sinh(x)ẋ = 0.

Let MARKER 9 represent the center of mass of the bead. Assume a penalty of 1E4 is used to enforce the constraint, and a penalty of 1E2 to enforce the time derivative of the constraint. The MotionSolve Force_Penalty statement is:

<Force_Penalty

id = "1"

label = "Particle sliding on a catenary"

type = "EXPRESSION"

expr = "dy(9) - cosh(dx(9))"

penalty = "1E4"

penalty1 ="1E2"

The displacement constraint for the penalty force is:

g = dy(9)-cosh(dx(9))

The velocity constraint for the penalty force is:

ġ = vy(9)-sinh(dx(9)*vx(9)

The penalty force that is computed by MotionSolve is:

F = -1E4*g - 1E2*ġ

= -1E4*[dy(9)-cosh(dx(9)] -1E2*[vy(9)-sinh(dx(9))*vx(9)]

The generalized force that is computed by MotionSolve along global-X and global-Y is: force_penalty_equ8

Example-2:

This example demonstrates how a simple unilateral constraint can be implemented with Force_Penalty. Figure 3 below shows a line L, y=1+x, defined in the coordinate system of a Marker, J.

force_penalty_fig3

Figure 3: A particle P bouncing on a line L: y=x+1

J has origin O, x-axis along X and y-axis along Y. A point P, belonging to another body, has coordinates (xp,yp) as measured in J. For the sake of convenience, define a Marker I=10, whose origin is at P.

The constraint, which we wish to impose, is that as the point P moves, it must not penetrate the material under L. The objectives of this example are:

•	Define the inequality constraint that captures this requirement.

•	Define a Force_Penalty to apply a contact force dependent on the penetration and penetration velocity.

From simple coordinate geometry considerations, we can see that:

•	The green region, of permissible motion, is is represented by: yp > 1 + xp

•	The red region, of impermissible motion, is represented by: yp < 1 + xp

•	The boundary between these regions is represented by line L: yp = 1 + xp

The constraint we wish to impose is: yp ≥ 1 + xp or 1 + xp – yp ≤ 0

The Force_Penalty definition is therefore:

<Force_Penalty

id = "10"

label = "Enforce 1+x < = y"

type = "EXPRESSION"

expr = "1 + dx(10) - dy(10)"

penalty = "1E4"

penalty = "1E2"

unilateral = "True"

smoothing_factor = "1.0"

Example-3:

This example demonstrates how a simple unilateral constraint can be implemented with Force_Penalty. Figure 3 below shows a line L, y=1+x, defined in the coordinate system of a Marker, J.

force_penalty_fig4

Figure 4: A particle P bouncing between two lines L1: z=1+x and L2: z=1-x

J has origin O, x-axis along X and Z-axis along Z. A point P, belonging to another body, has coordinates (xp,zp) as measured in J. For the sake of convenience, define a Marker I=10, whose origin is at P. The constraint, which we wish to impose, is that as the point P moves, it must not penetrate the material under either L1 or L2. The objectives of this example are:

•	Define the inequality constraints that captures this requirement.

•	Define two Force_Penalty to enforce these constraints.

From simple coordinate geometry considerations, we can see that:

•	The green region, of permissible motion above L1, is represented by: zp > 1 + xp

•	The red region, of impermissible motion below L1, is represented by: zp < 1 + xp

•	The boundary between these regions is represented by line L1: zp = 1 + xp

Similarly, we deduce that:

•	The green region, of permissible motion above L2, is represented by: zp > 1 - xp

•	The blue region, of impermissible motion below L2, is represented by: zp < 1 - xp

•	The boundary between these regions is represented by line L2: zp = 1 - xp

The constraints we wish to impose are:

•	zp ≥ 1 + xp or 1 + xp – zp ≤ 0

•	zp ≥ 1 - xp or 1 - xp – zp ≤ 0

The Force_Penalty definitions are therefore:

<Force_Penalty |<Force_Penalty

id = "1" | id = "2"

label = "Enforce 1+x<=z" | label = "Enforce 1-x<=z"

type = "EXPRESSION" | type = "EXPRESSION"

expr = "1+dx(10)-dz(10)" | expr = "1-dx(10)-dz(10)"

penalty = "1E4" | penalty = "1E4"

penalty1 = "1E2" | penalty1 = "1E2"

unilateral = "True" | unilateral = "True"

smoothing_factor = "1.0" | smoothing_factor = "1.0"

/> |/>

Figure 5 below shows the X and Z coordinates of particle P (Marker 10), when it starts from an initial condition of Xp=100 mm and Zp=120 mm and falls under the influence of gravity. The acceleration due to gravity is -9810 mm/s2 in the Z direction.

The particle bounces between the two lines and eventually comes to rest at Zp=1 and Xp=0.

force_penalty_fig5

Figure 5: Z and X time history of particle P bouncing between two lines L1: z=1+x and L2: z=1-x

Model Element
Description
PFORCE defines a penalty force whose purpose is to maintain a “soft” constraint in the system. This means that when the constraint is satisfied, there is no penalty force. However, if the constraint is violated, a penalty force that tries to reduce the violation is generated. The generalized force acts on all the coordinates involved in the definition of the constraint.
Declaration
def PFORCE(ID, LABEL="", FUNCTION="", PENALTY=0.0, PENALTY1=0.0, UNILATERAL=False, SMOOTHING_FACTOR=1.0, ROUTINE="", INTERPRETER="", SCRIPT=""):
Attributes
id	Element identification number (integer>0). This number is unique among all the PFORCE elements.
LABEL	This attribute describes the name of the PFORCE element. The description is primarily used to make the input more readable.
FUNCTION	Specifies the MotionSolve expression that defines the algebraic constraint. Any valid run-time MotionSolve expression can be provided as input. Or The list of parameters that are passed from the data file to the user defined subroutine, PFOSUB.
PENALTY	Specifies a penalty factor to be used in calculating the restoring force that is used to enforce that the algebraic constraint is always zero. Penalty ≥ 0.
PENALTY1	Specifies a second penalty factor to be used in calculating the restoring force that is used to enforce that the time derivative of the algebraic constraint. The default value for Penalty1 is zero. When specified, Penalty1 ≥ 0.
UNILATERAL	A Boolean attribute. Select one from TRUE and FALSE. TRUE indicates an inequality constraint (unilateral), and FALSE indicates an equality constraint (bilateral) is being specified. Unilateral is an optional attribute. It defaults to FALSE when not specified.
SMOOTHING_FACTOR	The penalty force due to a unilateral constraint is ramped up in a smooth manner using a STEP function. This specifies the x-value at which the smoothing is completed. SMOOTHING_FACTOR > 0.
ROUTINE	Specifies an alternative name for the user subroutine.
INTERPRETER	Specifies the interpreted language that the user script is written in. Valid choices are MATLAB or PYTHON.
SCRIPT	Specifies the path and name of the user written script that contains the routine.
Comments See Force_Penalty
Example The first example shows how a PFORCE element may be defined. PFORCE(1,LABEL="Particle sliding on a catenary", FUNCTION="DY(9)-COSH(DX(9))", PENALTY=1e4, PENALTY1=1e2) The second example shows how a unilateral PFORCE element may be defined. PFORCE(10,LABEL="Enforce 1+x < = y", FUNCTION="1+DX(10)-DY(10)" PENALTY=1e4, PENALTY1=1e2,UNILATERAL=True, SMOOTHING_FACTOR=1)

Force: Penalty

Force: Penalty

Force: Penalty

Model Element

Description

Format

Attributes

Comments

Examples

Model Element

Description

Declaration

Attributes

id

Element identification number (integer>0). This number is unique among all the PFORCE elements.

LABEL

This attribute describes the name of the PFORCE element. The description is primarily used to make the input more readable.

FUNCTION

PENALTY

Specifies a penalty factor to be used in calculating the restoring force that is used to enforce that the algebraic constraint is always zero. Penalty ≥ 0.

PENALTY1

Specifies a second penalty factor to be used in calculating the restoring force that is used to enforce that the time derivative of the algebraic constraint. The default value for Penalty1 is zero. When specified, Penalty1 ≥ 0.

UNILATERAL

A Boolean attribute. Select one from TRUE and FALSE. TRUE indicates an inequality constraint (unilateral), and FALSE indicates an equality constraint (bilateral) is being specified. Unilateral is an optional attribute. It defaults to FALSE when not specified.

SMOOTHING_FACTOR

ROUTINE

INTERPRETER

Specifies the interpreted language that the user script is written in. Valid choices are MATLAB or PYTHON.

SCRIPT

Specifies the path and name of the user written script that contains the routine.

Comments

See Force_Penalty

Example

The first example shows how a PFORCE element may be defined.

The second example shows how a unilateral PFORCE element may be defined.

See Also: