10罚函数法
第十章 罚函数法与广义乘子法本章主要内容：外罚函数法 内罚函数法 广义乘子法教学目的及要求：了解外罚函数法、内罚函数法、广义乘子法。
教学重点：外罚函数法． 教学难点：广义乘子法． 教学方法：启发式．教学手段：多媒体演示、演讲与板书相结合． 教学时间：3学时． 教学内容：§10.1 外罚函数法考虑约束非线性最优化问题min ();()0,1,2,,,()0,1,2,,,i j f x g x i m h x j l ⎧⎪≥=⎨⎪==⎩s.t. （10.1.1）其中(),()1,2,,()1,2,,i j f x g x i m h x j l ==（）和（）都是定义在n R 上的实值连续函数．记问题（10.1.1）的可行域为S ．对于等式约束问题min ();()0,1,2,,,j f x h x j l ⎧⎨==⎩s.t. （10.1.2）定义罚函数211(,)()()lj j F x f x h x σσ==+∑，其中参数σ是很大的正常数．这样，可将问题（10.1.2）转化为无约束问题 1min (,)nx RF x σ∈． （10.1.3） 对于不等式约束问题min ();()0,1,2,,,i f x g x i m ⎧⎨≥=⎩s.t. （10.1.4）定义罚函数221(,)()[max{0,()}]mi i F x f x g x σσ==+-∑，其中σ是很大的正数．这样，可将问题（10.1.4）转化为无约束问题2min (,)nx RF x σ∈． （10.1.5） 对于一般的约束最优化问题（10.1.1），定义罚函数 (,)()()F x f x P x σσ=+， 其中σ是很大的正数，()P x 具有下列形式：11()[()][()]mli j i j P x g x h x ϕφ===+∑∑．ϕ和φ是满足下列条件的实值连续函数：()0,0()0,0y y y y ϕϕ=≥⎧⎨><⎩当时，当时, ()0,0()0,0y y y y φφ==⎧⎨>≠⎩当时，当时.函数ϕ和φ的典型取法为()[max{0,}]y y αϕ=-， ()||y y βφ=，其中1α≥和1β≥是给定的常数，通常取作1或2．这样，可将约束问题（10.1.1）转化为无约束问题min (,)()()nx RF x f x P x σσ∈=+， （10.1.6） 其中σ是很大的正数，()P x 是连续函数．算法10-1（外罚函数法）Step1 选取初始数据．给定初始点0x ，初始罚因子10σ>，放大系数1α>，允许误差0ε>，令1k =．Step2 求解无约束问题．以1k x -为初始点，求解无约束问题min (,)()()nx RF x f x P x σσ∈=+， （10.1.11） 设其最优解为k x ．Step3 检查是否满足终止准则．若()k k P x σε<，则迭代终止，k x 为约束问题（10.1.1）的近似最优解；否则，令1,:1k k k k σασ+==+，返回Step2．例1 用外罚函数法在平面23x y +=上求一点(,)P x y ，使得P 到定点(3,2)A 的距离近似最短，取310.5,2,10σαε-===．解 令22(,)(3)(2)f x y x y =-+-，问题可归为求解如下最优化问题22min (,)(3)(2);(,)230.f x y x y h x y x y ⎧=-+-⎨=+-=⎩s.t.引入罚函数 222(,,)(3)(2)(23)k k F x y x y x y σσ=-+-++-． 则原约束最优化问题相应的一系列无约束最优化问题为：min (,,)k F x y σ，其中110.5,2,1,2,k k k σσσ+===．解上述无约束问题，得31122(,),1515TT k k k k k k x y σσσσ⎛⎫++= ⎪++⎝⎭，同时216(,)(15)kk k k k P x y σσσ=+． 依次对1,2,k =用上述公式计算(,)T k k x y 和(,)k k k P x y σ，结果如下表所示．由迭代终止条件3216(,)10(15)k k k k k P x y σσσ-=<+可得原约束问题的近似最优解（保留4位有效数字）1212(,)(2.2002,0.4003)T T x y =．§10.2 内罚函数法对于不等式约束问题（10.1.4），其可行域S 的内部int {|,()0,1,2,,}n i S x x R g x i m =∈>=．为了保持迭代点始终含于int S ，我们引入另一种罚函数(,)()()G x r f x rB x =+，其中r 是很小的正数，()B x 是int S 上的非负实值连续函数，当点x 趋向可行域S 的边界时，()B x →∞．显然，罚函数(,)G x r 的作用是对企图脱离可行域的点给予惩罚，相当于在可行域的边界设置了障碍，不让迭代点穿越到可行域之外，因此也称为障碍函数． (,)G x r 的两种常用的形式为11(,)()()mi i G x r f x r g x ==+∑， 及1(,)()ln ()mi i G x r f x r g x ==-∑分别称为倒数障碍函数和对数障碍函数．算法10-2（内罚函数法或SUMT 内点法）Step1 选取初始数据．给定初始内点0int x S ∈，初始参数10r >，缩小系数(0,1)β∈，允许误差0ε>，令1k =．Step2 求解无约束问题．以1k x -为初始点，求解无约束问题min (,)()();int ,k k G x r f x r B x x S =+⎧⎨∈⎩s.t. （10.2.2）设其最优解为k x ．Step3 检查是否满足终止准则．若()k k r B x ε<，则迭代终止，k x 为约束问题（10.1.4）的近似最优解；否则，令1,:1k k r r k k β+==+，返回Step2．算法10-3（求内罚函数法中初始内点的算法）Step1 选取初始数据．给定初始点0x ，初始参数10r >，缩小系数(0,1)β∈，令1k =．Step2 确定指标集．令{|()0,1,2,,}k i k I i g x i m =≤=，{|()0,1,2,,}k i k J i g x i m =>=．Step3 检验是否满足终止准则．若k I =∅，则int k x S ∈，迭代终止；否则，转Step4．Step4 求解无约束问题．以k x 为初始点，求无约束问题1min ()kk x S H x +∈的最优解1k x +，其中111()()()kkk i k i I i J i H x g x r g x ++∈∈=-+∑∑， {|()0,}k i k S x g x i J =>∈．令21,:1k k r r k k β++==+，返回Step2．§10.3 广义乘子法考虑只带等式约束的最优化问题（10.1.2）．不难知道，问题（10.1.2）等价于问题21min ()();2()0,1,2,,,l j j j f x h x h x j l σ=⎧⎡⎤+⎪⎢⎥⎨⎣⎦⎪==⎩∑s.t. （10.3.1）其中0σ>．该问题的Lagrange 函数211(,,)()()()2l ljj jj j L x v f x h x v h x σσ===+-∑∑， 由于(,,)L x v σ中既有罚项21()2ljj h x σ=∑，又有乘子项1()lj jj v h x =-∑，故称为乘子罚函数．与外罚函数类似，若设{}k σ为严格单调递增且趋于无穷大的正数列，则我们可以把求解等式约束问题（10.1.2）转化为求解一系列的无约束问题2()11min (,,)()()()2nl lk kk k j jj x Rj j L x v f x h x vh x σσ∈===+-∑∑， （10.3.2）其中()()()12(,,,)k k k T k l v v v v =是第k 次迭代中采用的Lagrange 乘子，并且有定理10.3.1 设k x 是无约束问题（10.3.2）的最优解，则k x 也是问题min ();()(),1,2,,j j k f x h x h x j l⎧⎨==⎩s.t. （10.3.3）的最优解．记12()((),(),,())T l h x h x h x h x =．算法10-4（等式约束下的广义乘子法）Step1 选取初始数据．给定初始点0x ，初始乘子1v ，初始罚因子10σ>，放大系数1α>，允许误差0ε>，参数(0,1)γ∈，令1k =．Step2 求解无约束问题．以1k x -为初始点，求解无约束问题（10.3.2），设其最优解为k x ．Step3 检查是否满足终止准则．若()k h x ε<，则迭代终止，k x 为等式约束问题（10.1.2）的近似最优解；否则，转Step4．Step4 判断收敛快慢．若1()()k k h x h x γ-≥，则令1k k σασ+=，转Step5；否则，令1:k k σσ+=，转Step5．Step5 进行乘子迭代．令(1)()(),1,2,,,:1k k j j k j k v v h x j l k k σ+=-==+，返回Step2．。
物流配送中几种路径优化算法
捕食搜索算法 动物学家在研究动物的捕食行为时发现，尽管由于动物物种的不同而造成的身体结构的千差万别，但它们的捕食行为却惊人地相似.动物捕食时，在没有发现猎物和猎物的迹象时在整个捕食空间沿着一定的方向以很快的速度寻找猎物.一旦发现猎物或者发现有猎物的迹象，它们就放慢步伐，在发现猎物或者有猎物迹象的附近区域进行集中的区域搜索，以找到史多的猎物.在搜寻一段时间没有找到猎物后，捕食动物将放弃这种集中的区域，而继续在整个捕食空间寻找猎物。
模拟动物的这种捕食策略，Alexandre于1998提出了一种新的仿生计算方法，即捕食搜索算法(predatory search algorithm, PSA)。
基本思想如下：捕食搜索寻优时，先在整个搜索空间进行全局搜索，直到找到一个较优解;然后在较优解附近的区域(邻域)进行集中搜索，直到搜索很多次也没有找到史优解，从而放弃局域搜索;然后再在整个搜索空间进行全局搜索.如此循环，直到找到最优解(或近似最优解)为止，捕食搜索这种策略很好地协调了局部搜索和全局搜索之间的转换.目前该算法己成功应用于组合优化领域的旅行商问题(traveling salesm an problem )和超大规模集成电路设计问题(very large scale integrated layout)。
捕食搜索算法设计 (1)解的表达 采用顺序编码，将无向图中的，n一1个配送中心和n个顾客一起进行编码.例如，3个配送中心，10个顾客，则编码可为:1一2一3一4一0一5一6一7一0一8一9一10其中0表示配送中心，上述编码表示配送中心1负贡顾客1,2,3,4的配送，配送中心2负贡顾客5,6,7的配送，配送中心3负贡顾客8,9,10的配送.然后对于每个配送中心根据顾客编码中的顺序进行车辆的分配，这里主要考虑车辆的容量约束。
依此编码方案，随机产生初始解。
(2)邻域定义 4 仿真结果与比较分析(Simulation results and comparison analysis) 设某B2C电子商务企业在某时段由3个配送中心为17个顾客配送3类商品，配送网络如图2所示。
第二节 罚函数法
β ∈ (0,1) ，允许误差 ε > 0 ，置 k = 1 ；
k −1 x 为初点，求解无约束优化问题 Step2: 以
min G ( x, r ) = f ( x) + rB( x) s.t. x ∈ int S
设其极小点为 x ；
G ( x, r ) = f ( x) − r ∑ ln g i ( x)
m
-----对数障碍函数
由 G( x, r ) 的定义， r 取值越小，问题
min G ( x, r ) = f ( x) + rB( x) s.t. x ∈ int S
的最优解越接近约束优化问题的最优解。 2. 内点罚函数法的计算步骤
k min H ( x) x Step4: 以 为初始点求解无约束问题 x∈Sk k +1 的最优解
x k +1 ，其中
H k +1 ( x) = −∑ gi ( x) + rk +1 ∑ gi ( x)
i∈I k i∈J k
， Sk = {x | gi ( x) > 0, i ∈ J k }
令 rk + 2 = βrk +1 , k = k + 1, 返回 Step2. 注：该算法中，对于 k = 0,1,2,L ，有 I k +1 ⊂ I k , J k +1 ⊃ J k ，且 最后某个 I k = ∅ 。 三. 广义乘子法 1. 对于等式约束优化问题
φ ( x, y , ω , σ ) = f ( x ) − ∑ ω j ( g j ( x ) − y j ) +
2 j =1 l
罚函数法(SUMT法)(ppt文档)
( NP) 求解 min( X , M ) 设其最优解为 X*(M)， XRn
研究 X*(M) 与(NP)的最优解 X* 之间的关系
线性规划3-6
一.外点法迭代原理
(NP) min f (X )
m
( X , M ) f ( X ) M [min( 0, gi ( X ))]2
i 1
gi (X ) 0 gi (X )
30 若 gi ( X (k ) ( Mk )) , i 1,2, ,m, 则迭代终止，X X (k ) ( M k ) 否则取Mk+1=C Mk , 其中C = 5~10
D
X
X(k)(Mk )
gi (X ) 0
令 k:= k+1 转20
M4 1000 X (4) (1000)
D
g2(X ) 0
g1( X ) 0
X(k)(Mk )
X
x1
线性规划3-6
一.外点法迭代原理
(NP) min f (X )
s.t. gi ( X ) 0
线性规划3-6
第三章 非线性规划
一.外点罚函数法(外点法)
( NP) min f (X )
s.t
.
gi hj
( (
X X
) )
0, 0,
i j
1, 2, 1, 2,
,m ,p
外点法迭代原理 外点法迭代步骤 外点法举例 外点法的优缺点
二.外点法迭代步骤
(NP) min f (X )
通过迭代逐渐增大罚因子M：
(NP)min f (X )
s.t. gi ( X ) 0
《惩罚函数法》PPT课件
这样的存在迫使最优解在可行域内取得。 随着的增大或更特殊地取为+∞，则问题(*)就成为：
min (x12+x22) 当(x1+x2-2)=0.
这恰为所要求解的原问题.
引例求解思想的理论支持
问题 min (x12+x22)+(x1+x2-2)2
最优解的解析式为：
x1() x2 () 221
一般地，对于等式约束问题， min f(x) s.t. hj(x)=0, j=1:n
将此问题改造成一个新问题(**)：
n
min F( x),f(x )h2 j(x),其 中 为一个大正数 j1
这个新问题的最优解 必定使～x得hj( )接近于0 ～x 否则的话式子中的第二项就会是一个很大的正数 现在的这个点 就不～会x 是这个无约束问题的极小点
[max2 {(x0 )},2-]s.....[.maxm {(x 0),}2 -]s
h12(x)h22(x).....h .n2(x)）
m
n
f(x )( [mai( x x ){ } 2 0 ] ,h -j2 s (x ))
i 1
j 1
P( x)
F(x, )-----增广目标函数
P(x)-----惩罚函数（惩罚项） ----罚因子
min x12+x22 s.t. x1+x2-2=0
由图解法易见最优解为(1,1)T
将这个问题改造为一个无约束问题如下：
min (x12+x22)+(x1+x2-2)2 (*)
为一个充分大的正的参数
min x12+x22 s.t. x1+x2-2=0
最优化方法 第三章(罚函数法)
这种惩罚策略，对于在无约束的求解过程中企图违反约
束的迭代点给予很大的目标函数值，迫使无约束问题的 极小点或者无限地向可行域D靠近，或者一直保持在可 行域D内移动，直到收敛到原来约束最优化问题的极小 点。
不改变可行域局部极小值，可以将 约束域之外的局部极小值变大。
p ( x) 0, x D p ( x) 0, x D
k k
k 1
k 1
xk 1是F x, M k 1 的最优解.
k 1 k k 1 k 0 M k 1 M k p ( x ) p ( x ) p ( x ) p ( x )
M k 1 M k
(3) f ( x k 1 ) M k p( x k 1 ) F ( x k 1 , M k ) F ( x k , M k ) f ( x k ) M k p( x k )
gi ( x) gi ( x) max gi ( x), 0 = 罚函数p(x)的构造 2 m l p( x) (max gi ( x), 0) 2 h 2 j ( x)
i 1 j 1
(1) p(x)连续 (2) p( x) 0, x D (3) p( x) 0, x D
二、外点法 外点罚函数法算法步骤 1：给定初始点 x 0 ，初始罚因子M1 0 (可取M1 1 ), 精度 0, k : 1. 2：以 x k 1初始点，求解无约束优化问题
min F ( x, M k ) f ( x) M k p( x)
得到极小点 x* ( M k )，记为 x k ， 其中
p( x) (max gi ( x), 0) h 2 j ( x)
2 i 1 j 1 m l
非线性规划问题的求解方法
输入参数语法：
x = fmincon(fun,x0,A,b) x = fmincon(fun,x0,A,b,Aeq,beq) x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub) x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon) x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options) x= fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options,P1,P2, ...)
4、其它求解算法
（1）间接法 （2）直接法
直接搜索法 以梯度法为基础的间接法
无约束规划的Matlab求解函数 数学建模案例分析（截断切割，飞机排队）
（1）间接法
在非线性最优化问题当中，如果目标函 数能以解析函数表示，可行域由不等式约束 确定，则可以利用目标函数和可行域的已知 性质，在理论上推导出目标函数为最优值的 必要条件，这种方法就称为间接法（也称为
第二步：求 (k) 最优的目标函数
function r=fungetlamada(lamada) %关于lamada的一元函数，求最优步长 global x0 d=fun1gra(x0); r=2*(x0(1)-lamada*d(1))^2+(x0(2)lamada*d(2))^2; %注意负号表示是负梯度
三、Matlab求解有约束非线性规划
1. 用fmincon函数求解形如下面的有约束 非线性规划模型
一般形式：
min f ( X ) s.t. AX b
Aeq X beq l X u c(X ) 0 ceq ( X ) 0
最优化方法之-罚函数法讲解
内点罚函数法优点
迭代总在可行域内进行，每一个中间结果都是 可行解，可以作为近似解。
内点罚函数法缺点
选取初始可行点较困难，且只适用于含不等式 约束的非性性规划问题。
x2)
0
x1
1
2 2M
x2
1
(2 2M )2
1
2M
当M
时，有
x1 x 2
0 0
步骤:
1.给定x( 初 0) ，始 初点 始 M10 罚 (M1 因 1 ) ， 子 放 系c数 1 ，允许 0 ， 误 k置 差 1 。
2.以x(k1)为初始点，求解无约 问题 束 minf(x) Mkp(x)
201 101
以此类推，得序列： 3，21，201，2001， 2 11 101 1001
2
外点罚函数法的一个重要特点:
函数 F(x,M)是在整个En空 内间 进行,优 初化 始
点可任意,且 选外 择点法也可规 用划 于的 非最 凸.优
缺点:
1.惩罚项Mp(x)的二阶偏导数一般在 不;存
2.外点法的中间可 结行 果,不 解 不能 是作为近 ;
P~(x,rk) iSk
gi(x)rk
iTk
1 gi(x)
记 R~k x|gi(x)0 iTk
(rk 0)
5.以x(k)为初始点R， ~k域在内，求障碍 P~(函 x,rk数 ) 的极小:点
minP~(x,rk) s.t. xR~k 得x(k1)，转 6。
6.令 0rk1rk如r取 k1110rk,置 k:k1， 转 2。
引入罚项
p(x )
m
g i( x )
i 1
其中 ( y )是连续函数，且满足
(y ) 0 ( y ) 0
函数极限的十种求法
函数极限的十种求法信科2班江星雨250函数极限可以分成而运用ε-δ定义更多的见诸于已知极限值的证明题中。
掌握这类证明对初学者深刻理解运用极限定义大有裨益。
以的极限为例，f(x) 在点以A为极限的定义是：对于任意给定的正数ε（无论它多么小），总存在正数，使得当x满足不等式时，对应的f(x)函数值都满足不等式：，那么常数A就叫做函数f(x)当x→x。
时的极限。
1.利用极限的四则运算法则：极限四则运算法则的条件是充分而非必要的，因此，利用极限四则运算法则求函数极限时，必须对所给的函数逐一进行验证它是否满足极限四则运算法则条件，满足条件者。
方能利用极限四则运算法则进行求之。
不满足条件者，不能直接利用极限四则运算法则求之。
但是，井非不满足极限四则运算法则条件的函数就没有极限，而是需将函数进行恒等变形，使其符合条件后，再利用极限四则运算法则求之。
而对函数进行恒等变形时，通常运用一些技巧如拆项、分子分母同时约去零因子、分子分母有理化、通分、变量替换等等。
例 1求lim( x 2 − 3x + 5).x→ 2解：lim( x 2 − 3x + 5) = lim x 2 − lim 3x + lim 5= (lim x) 2 − 3 lim x + lim 5= 2 2 − 3 ⋅ 2 + 5 = 3.x→2 x →2 x →2 x →2 x →2 x →2 x →22.利用洛必达法则洛必达(L 'Hopital)法则是在一定条件下通过分子分母分别求导再求极限来确定未定式值的方法.简单讲就是，在求一个含分式的函数的极限时，分别对分子和分母求导，在求极限，和原函数的极限是一样的。
一般用在求导后为零比零或无穷比无穷的类型。
利用洛必达求极限应注意以下几点：设函数f(x)和F(x)满足下列条件：（1）x→a时,lim f(x)=0,lim F(x)=0;（2）在点a的某去心邻域内f(x)与F(x)都可导,且F(x)的导数不等于0;（3）x→a时,lim(f'(x)/F'(x))存在或为无穷大则x→a时,lim(f(x)/F(x))=lim(f'(x)/F'(x))例1：1-cosx = 1-{1-2[sin(x/2)]^2} = 2[sin(x/2)]^2xsinx = 2xsin(x/2)cos(x/2)原式= lim 2[sin(x/2)]^2 / [2xsin(x/2)cos(x/2)] = tgx / x对分子分母同时求导（洛必达法则）(tgx)' = 1 / (cosx)^2(x)' = 1原式= lim 1/(cosx)^2当x --> 0 时,cosx ---> 1原式= 13.利用两个重要极限：应用第一重要极限时，必须同时满足两个条件：①分子、分母为无穷小，即极限为0 ；②分子上取正弦的角必须与分母一样。
接触问题的罚函数法
BULETINUL INSTITUTULUI POLITEHNIC DIN IAŞIPublicat deUniversitatea Tehnică …Gheorghe Asachi” din IaşiTomul LIV (LVIII), Fasc. 3, 2011SecţiaCONSTRUCŢII. ARHITECTURĂPENALTY BASED ALGORITHMS FOR FRICTIONALCONTACT PROBLEMSBYANDREI-IONUŢŞTEFANCU*, SILVIU-CRISTIAN MELENCIUCand MIHAI BUDESCU“Gheorghe Asachi” Technical University of Iaşi,Faculty of Civil Engineering and Building ServicesReceived: June 21, 2011Accepted for publication: August 29, 2011Abstract. The finite element method is a numerical method that can be successfully used to generate solutions for problems belonging to a vast array of engineering fields: stationary, transitory, linear or nonlinear problems. For the linear case, computing the solution to the given problem is a straightforward process, the displacements are obtained in a single step and all the other quantities are evaluated afterwards. When faced with a nonlinear problems, in this case with a contact nonlinearity, one needs to account for the fact that the stiffness matrix of the systems varies with the loading, the force vs. stiffness relation being unknown prior to the beginning of the analysis. Modern software using the finite element method to solve contact problems usually approaches such problems via two basic theories that, although different in their approaches, lead to the desired solutions. One of the theories is known as the penalty function method, and the other as the Lagrange multipliers method. The hereby paper briefly presents the two methods emphasizing the penalty based ones. The paper also underscores the influence of input parameters for the case of the two methods on the results when using the software ANSYS 12.Key words: finite element method; pure penalty methods; Lagrange multipliers method; H-adaptive meshing.*Corresponding author: e-mail: stefancu@ce.tuiasi.ro120Andrei-IonuţŞtefancu, Silviu-Cristian Melenciuc and Mihai Budescu1. IntroductionThe finite element method (FEM) is a numerical method that can be applied to obtain solutions to problems belonging to a variety of engineering disciplines: stationary problems, transitory problems, linear or nonlinear stress/strain problems. Heat transfer, fluid flow, or electromagnetism problems can also be solved using FEM.The first ones to publish articles in this field are Alexander Hrennikoff (1941) and Richard Courant (1943). Although their approaches are different, they have a common feature – discretization of continua into a series of discrete sub-domains called elements. Olgierd Zienkiewicz summarizes the work of his predecessors in what was to be known as FEM (1947). In 1960 Ray Clough is the first one to use the term finite element. Zienkiewicz and Cheung published in 1967 the first book entirely devoted to the finite element method.In this context, in 1971, the first version of the ANSYS software is released. Currently, the ANSYS software package is able to solve static, dynamic, heat transfer, fluid flow, electromagnetics problems, etc. ANSYS is a market leader for more than 20 years (Moavenim, 1999). In the present study version 12.0.1 of the software was used (ANSYS Workbench 2.0, 2009).2. Finite Element Formulations of Contact ProblemsThere are two basic theories that, although different in their approaches, offer the desired solutions to body contact problems: the penalty function method and the Lagrange multipliers method.The main difference between them is the way they include in their formulation the potential energy of contacting surfaces.The penalty function method, due to its economy, has received a wider acceptance. The method is very useful when solving frictional contact problems, while the Lagrange method, based on multipliers, is known for its accuracy.The main drawback of the Lagrange method is that it may lead to ill-converging solutions while the penalty formulation may lead to inaccurate ones.In the following the pure penalty, the augmented Lagrange methods will be presented.2.1. The Penalty MethodThe penalty method involves adding a penalty term to enhance the solving process. In contact problems the penalty term includes the stiffness matrix of the contact surface. The matrix results from the concept that one body imaginary penetrates the another (Wriggers et al., 1990).Bul. Inst. Polit. Ia şi, t. LVII (LXI), f. 3, 2011 121The stiffness matrix of the contact surface is added to the stiffness matrix of the contacting body, so that the incremental equation of the Finite Element becomes[K b + K c ]u = F , (1)where: K b is the stiffness matrix of contacting bodies; K c – stiffness matrix of contact surface; u – displacement; F – force.The magnitude of the contact surface is unknown (Stein & Ramm, 2003), therefore its stiffness matrix, K c , is a nonlinear term. The total load and displacement values are∑Δ=F F tot , (2)∑Δ=u u tot, (3)where: F tot is the force vector; u tot – displacement vector.To derive the stiffness matrix, the contact zone (encompassing the contact surface) is divided into a series of contact elements. The element represents the interaction between the surface node of one body with the respective element face of the other body. Fig. 1 shows a contact element in a two dimensional application. It is composed of a slave node (point S ) and a master line, connecting nodes 1 and 2. S 0 marks the slave node before the application of the load increment, and S marks the node after loading.Fig. 1 – Contact element – penalty method formulation.Given the nature of the numerical simulations presented afterwards only the sliding mode of friction will be presented. In this case, the tangential force acting at the contact surface equals the magnitude of the friction force, hence the first variation of the potential energy of a contact element ist n n d t t n n t t n n c g g k g g g k g f g f δμδδδδ)sgn( +=+=Π, (4)122 Andrei-Ionu ţ Ştefancu, Silviu-Cristian Melenciuc and Mihai Budescuwhere: k n represents penalty terms used to express the relationship between the contact force and the penetrations along the normal direction; k t – penalty terms used to express the relationship between the contact force and the penetrations along the tangential direction; g n – penetration along the normal direction; g t – penetration along the tangential direction;n n n f k g =, (5))()sgn(n n d t t g k g f μ−=. (6)2.2. The Augmented Lagrange Multiplier MethodIn the case of classical Lagrange Multiplier Method the contact forces are expressed by Lagrange multipliers. The augmented Lagrange method involves the regularization of classical Lagrange method by adding a penalty function from the penalty method (Simo & Laursen, 1992). This method, unlike the classical one, can be applied to sticking friction, sliding friction, and to africtionless contactFig. 2 – Contact element – Lagrangian methods.The contact problem involves the minimization of potential П by equating to zero the following expression:kg g g u u T T b 21)(),(+Λ+Π=ΛΠ, (7) where ⎥⎥⎦⎤⎢⎢⎣⎡⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧=Λk t k n t n t n Tλλλλλλ,...,,2211, (8) with:n λ – Lagrange multiplier for the normal direction;t λ– Lagrange multiplier for the tangential direction;Bul. Inst. Polit. Ia şi, t. LVII (LXI), f. 3, 2011 12312k n n n 12k t t t g g g g ,,...,g g g ⎡⎤⎧⎫⎧⎫⎧⎫⎪⎪⎪⎪⎪⎪=⎢⎥⎨⎬⎨⎬⎨⎬⎢⎥⎪⎪⎪⎪⎪⎪⎩⎭⎩⎭⎩⎭⎣⎦. (9)3. Parametric Analysis of Frictional ContactIn order to illustrate the way that the contact algorithms may influence the results a parametric analysis is performed. The purpose of this analysis is to exemplify how various input parameters can alter the results.3.1. Finite Elements Formulations of Contact Problems in ANSYSThe Finite Elements (FE) software ANSYS, for the penalty method, assumes that contact force along the normal direction is written as follows:cont.penetr.cont.F x K Δ=Δ, (10)where: K cont. is the contact stiffness, defined by real constant FKN for the 17x contact elements (in the current analysis the 174 contact element is used); x penetr. – distance between two existing nodes on separate contact bodies; F cont. – contact force. ANSYS automatically chooses the real constant FKN as a scale factor of the stiffness of the underlying elements. This value can be modified by the user (via FKN – a scale factor).Given the fact that the augmented Lagrange method is actually a penalty method with penetration control, the contact force is computed according to eq. (10), the only difference being the contact stiffness formulationpenetr.cont.1x K i i +=+λλ, (11)where λi is a Lagrange multiplier.Although the Lagrange multipliers are condensed out at the element level, one can think regarding this method as the same as a regular penalty one except that the contact stiffness is “updated” per contact element (Imaoka, 2001).Similar to the normal direction, a real constant – FKT models the tangential stiffness of the contact.3.2. Adaptive Solutions in ANSYSIn order to overcome the influence of the meshing upon the final results of the analysis and to improve the accuracy of the solution an adaptive solution will be used.124 Andrei-Ionu ţ Ştefancu, Silviu-Cristian Melenciuc and Mihai BudescuIn ANSYS the desired accuracy of a solution can be achieved by means of adaptive and iterative analysis, whereby h -adaptive methodology is employed.The h -adaptive method begins with an initial FE model that is refined over various iterations by replacing coarse elements with finer ones in selected regions of the model. This is effectively a selective remeshing procedure.The criterion for which elements are selected for adaptive refinement depends on geometry and, for the current analysis, on a 10% allowable difference between the maximum values of the frictional (obtained in two consecutive runs with different meshes).The user-specified accuracy is achieved when convergence is satisfied as follows:)in ...,3,2,1( ,1001R n i E i i i −=<⎟⎠⎞⎜⎝⎛−+φφφ, (11) where:φis the result quantity; E – expected accuracy (10% for this case); R – the region on the geometry that is being subjected to adaptive analysis (entire geometry in this case); i – the iteration number.The results are compared from iteration i to iteration i + 1. Iteration in this context includes a full analysis in which h -adaptive meshing and solving are performed.For this case of adaptive procedures, the ANSYS product identifies the largest elements, which are deleted and replaced with a finer FE representation (ANSYS, 2009).The overall results show a good behavior of the model. Only two iteration are performed in order to satisfy reach the expected accuracy of the solution.Table 1 h-Adaptive Methodology Convergence HistoryIteration Frictional Stress, [MPa]Change, [%] Nodes Elements 1 0.130042,518 352 2 0.12739–2.0571 14,158 8,2343.3. The Model Used in the Parametric AnalysisThe model used, represented in Fig. 3, comprises two solids made up of nonlinear structural steel materials. The larger solid has its lower surface fixed while at the upper end interacts with the smaller solid via a frictional contact (coefficient of friction 0.2). A normal pressure of 0.5 MPa is applied on top of the smaller solid, and displacement is applied on the left hand side face.Bul. Inst. Polit. Iaşi, t. LVII (LXI), f. 3, 2011 125Fig. 3 – The model used in the parametric analysis. A – frictionalcontact; B – applied pressure; C – fixed support; D – applieddisplacement.Input parameters:a) FE formulation (P1); this parameter can take two values: 0 for augmented Lagrange method and 1 for pure penalty method;b) the normal contact stiffness factor FKN (P2) that varies between 0.01 and 1;c) the tangent contact stiffness factor FKT (P3) that varies between 0.01 and 1.Output parameters: a number of output parameters have been monitored, such as: maximum (P5) and minimum (P8) normal elastic strain, maximum (P9) and minimum (P10) shear elastic strain, maximum (P12) and minimum (P13) normal stress, maximum (P14) and minimum (P15) shear stress, maximum (P11) frictional stress, maximum (P6) penetration, analysis run time (P7), maximum stiffness energy (P16).3.4. Results of the Parametric AnalysisThe parametric analysis provides a wide range of information regarding the dependence of the output parameters on the input ones. Based on the relevance of the results only a limited amount of them will be presented.The local sensitivity chart allows one to appreciate the impact of the input parameters on the output ones. This means that the output is computed based on the change of each input independently of the current value of each input parameter. The larger the change of the output, the more significant is the input parameter that was varied (ANSYS, 2009). Since the local sensibilities can only be computed for continuous parameters (P1 is a discrete one) the sensibility chart will be presented for the pure penalty and augmented Lagrange method individually.126Andrei-IonuţŞtefancu, Silviu-Cristian Melenciuc and Mihai BudescuFig. 4 – Local sensibility chart – Lagrangian method.Fig. 5 – Local sensibility chart – pure penalty method.Bul. Inst. Polit. Iaşi, t. LVII (LXI), f. 3, 2011 127It can be seen, from Fig. 5, that the pure penalty method is less sensitiveto contact normal and tangent stiffness than the augmented Lagrange method.The only output parameter influenced by the contact stiffness is the analysisrun-time.In Figs. 4 and 5 only sensitivities of the three parameters (P5, P6, P7)have been presented because the sensitivities of the other are zero or almostzero. Based on this the variation of P5, P6 and P7, with P2 and P3, arepresented in what follows.a bFig. 6 – Variation of P5 with P2 and P3: a – augmented Lagrangemethod; b – pure penalty method.a bFig. 7 – Variation of P6 with P2 and P3: a – augmented Lagrangemethod; b – pure penalty method.As it can be seen from Figs. 6 and 7 there isn’t an exact pattern of thevariation of the maximum normal elastic strain (P5), maximum penetration128Andrei-IonuţŞtefancu, Silviu-Cristian Melenciuc and Mihai Budescu(P6) or analysis run time (P7) with the normal (P2) and tangent (P3) contactstiffness factor.Give the kinematic nature of the problem, and the contact type(frictional) it can be observed from Fig. 8 that the analysis run time (the timeneeded to compute a solution for the given problem) is tangent stiffnessdependent.a bFig. 8 – Variation of P7 with P2 and P3: a – augmented Lagrangemethod; b – pure penalty method.4. ConclusionsGiven the high nonlinear characteristic of the frictional contacts an extraattention is necessary to be paid to contact algorithms and their inputparameters. In such case an h-adaptive solution is recommended to be usedbecause such approach can “fade out” the influence of contact parameters onmost of the output parameters.If working circumstances require fulfilling certain limitations, accuracyconditions may be enforced, thus improving the confidence level of the finalsolution. One must keep in mind though, that an increased number of accuracyconvergence conditions leads to prohibitive analysis run time.REFERENCESClough R.W., The Finite Element Method in Plane Stress Analysis. Proc. 2nd ASCEConf. on Electron. Comp., Pittsburg,Pennsylvania, 1960.Courant R., Variational Methods for the Solution of Problems in Equilibrum and Vibra-tions. Bul. of the Amer. Mathem. Soc., 49,1-23 (1943).Hrennikoff A., Solution of Problems of Elasticity by the Frame-Work Method. ASME,J. of Appl. Mech., 8, 169–175 1(941).Bul. Inst. Polit. Ia şi, t. LVII (LXI), f. 3, 2011 129 Imaoka S., Sheldon’s ANSYS Tips and Tricks: Understanding Lagrange Multipliers .available on-line at /ansys/tips_sheldon/STI07_Lagrange_Mul-tipliers.pdf , 2001.Moavenim A., Finite Element Analysis – Theory and Applicaion with ANSYS . PrenticeHall, Upper Saddle River, New Jersey, 1999.Simo J.C., Laursen T.A., An Augmented Lagrangian Treatment of Contact ProblemsInvolving Friction . Comp. a. Struct., 42, 1, 97-116 (1992).Stein E., Ramm E., Error-Controlled Adaptive Finite Elements in Solid Mechanics . J.Wiley a. Sons, NY, 2003.Wriggers P., Vu Van T., Stein E., Finite Element Formulation of Large DeformationImpact-Contact Problems with Friction . Comp. a. Struct., 37, 3, 319-331 (1990).Zienkiewicz O.C., Cheung Y.K., The Finite Element Method in Continuum andStructural Mechanics . McGraw Hill, NY, 1967.Zienkiewicz O.C., The Stress Distribution in Gravity Dams . J. Inst. Civ. Engng., 27,244-271 (1947). * * * ANSYS Workbench 2.0 Framework Version: 12.0.1. ANSYS Inc., 2009. * * * Design Exploration . ANSYS Inc., 2009. * * * Theory Reference, Analysis Tools, ANSYS Workbench Product Adaptive Solutions . ANSYS Inc., 2009.ALGORITMI BAZA ŢI PE METODA COREC ŢIILOR UTILIZA ŢI ÎNREZOLVAREA PROBLEMELOR DE CONTACT CU FRECARE(Rezumat)Metoda elementului finit este o metod ă numeric ă ce poate fi aplicat ă cu succes pentru a ob ţine solu ţiile problemelor dintr-o multitudine de discipline inginere şti: probleme sta ţionare, probleme tranzitorii, liniare sau neliniare. În cazul liniar g ăsirea solu ţiei unei probleme date este un proces simplu. Deplas ările sunt ob ţinute într-un singur pas de analiz ă, tensiunile şi deforma ţiile fiind evaluate ulterior. În cazul problemelor neliniare – în acest caz neliniaritate de contact – trebuie s ă se ţin ă cont de faptul c ă matricea de rigiditate a sistemului variaz ă func ţie de înc ărcare, rela ţia for ţă vs . rigiditate nefiind cunoscut ă a priori . Programele moderne, ce folosesc metoda elementului finit pentru a rezolva probleme de contact, abordeaz ă de obicei astfel de probleme prin intermediul a dou ă teorii care, de şi diferite în abord ările lor, conduc la solu ţia dorit ă. Una dintre teorii este cunoscut ă sub numele de metoda corec ţiilor, iar cealalt ă ca metoda multiplicatorilor Lagrange. În lucrare se prezint ă pe scurt cele dou ă metode, accentul punându-se pe metodele bazate pe corec ţii. Lucrarea eviden ţiaz ă, de asemenea, influen ţa parametrilor de intrare caracteristici algoritmilor de rezolvare a problemelor de contact asupra rezultatelor atunci când se utilizeaz ă pachetul software ANSYS 12.。
惩罚函数法算法
x L(x, k , v k ) f (xk ) v (jk ) k h j (xk ) h j (xk ) 0
j 1
l
v(jk 1) v(jk ) k h j (xk ), j 1, 2,
,l
转化求解法（二）：增广乘子法
等式约束下的增广乘子法 Step1 选取初始数据。给定初始点 x 0 ，初始乘子 λ1 ，初始罚因子 1 0 ，放 大系数 1 ，允许误差 0 ，参数 (0, 1)，令 k 1 。
利用序列无约束极小化方法（SUMT）
min G(x, rk ) f (x) rk B(x) s.t. x int S
转化求解法（一）：罚函数法
内罚函数法
Step1 选取初始数据。给定初始点 x0 int S ，初始参数
(0,1)
允许误差
0 ，令 k 1
1 0 ，缩小系数
,m ,l
(1)
其中
f (x) ，gi (x)(i 1, 2,
, m) 和 h j (x)( j 1, 2,
, l)
都是定义在
min s.t.
R n 上的实值函数。记问题（1）的可行域为 S
,m ,l
。
f ( x) gi (x) 0, i 1, 2, h j (x) 0, j 1, 2,
基本是想
把罚函数与Lagrange函数结合起来，构造出更合适的新目标函数， 使得在罚因子适当大的情况下，借助于Lagrange乘子就能逐步达 到原约束问题的最优解。 由于这种方法要借助于Lagrange乘子的迭代进行求解 而又区别于经典的Lagrange乘子法，故称为广义乘子法。