化妆品网站模板下载,临沂做商城网站,蛋糕网站模版,建设网站需要哪些东西文章目录 非线性规划基本概念与结论凸集与凸函数极值条件无约束条件的极值判断条件有约束条件的极值判断条件 无约束非线性规划一维搜索算法步骤示例特点代码模板 最速下降法算法详细步骤 代码实现示例最优步长的求解 黄金分割法斐波那契法牛顿法阻尼牛顿法模式搜索法Powell方法… 文章目录 非线性规划基本概念与结论凸集与凸函数极值条件无约束条件的极值判断条件有约束条件的极值判断条件 无约束非线性规划一维搜索算法步骤示例特点代码模板 最速下降法算法详细步骤 代码实现示例最优步长的求解 黄金分割法斐波那契法牛顿法阻尼牛顿法模式搜索法Powell方法 有约束非线性规划外点罚函数法算法步骤Python 代码实现注意事项内点罚函数法原理倒数障碍函数对数障碍函数 算法步骤Python 代码示例使用对数障碍函数代码说明注意事项 非线性规划
基本概念与结论
有约束非线性规划的一般数学模型为如下形式 min f ( x ) s . t . { g i ( x ) ⩾ 0 , i 1 , 2 , . . . , m h j ( x ) 0 , j 1 , 2 , . . . , n \min f\left( x \right) \\ s.t.\left\{ \begin{array}{c} g_i\left( x \right) \geqslant 0,i1,2,...,m\\ h_j\left( x \right) 0,j1,2,...,n\\ \end{array} \right. minf(x)s.t.{gi(x)⩾0,i1,2,...,mhj(x)0,j1,2,...,n 下面给出非线性规划的基本概念 定义1设 f ( x ) f(x) f(x)为目标函数 S S S为可行域 x ∗ ∈ S x^*\in S x∗∈S。若 ∀ x ∈ S , f ( x ) ⩾ f ( x ∗ ) \forall x\in S,f\left( x \right) \geqslant f\left( x^* \right) ∀x∈S,f(x)⩾f(x∗)则称 x ∗ x^{*} x∗为 min f ( x ) s . t . x ∈ S \min f\left( x \right) \\ s.t. x\in S minf(x)s.t.x∈S 极小化问题的最优解(全局最优解)。 定义2设 f ( x ) f(x) f(x)为目标函数 S S S为可行域 x ∗ ∈ S x^*\in S x∗∈S。若 ∃ x ∗ ∈ { x ∣ ∥ x − x ∗ ∥ ε , ε 0 } , f ( x ) ⩾ f ( x ∗ ) \exist x^*\in \left\{ x\left| \left\| x-x^* \right\| \right. \varepsilon ,\varepsilon 0 \right\} ,f\left( x \right) \geqslant f\left( x^* \right) ∃x∗∈{x∣∥x−x∗∥ε,ε0},f(x)⩾f(x∗)则称 x ∗ x^{*} x∗为极小化问题 min f ( x ) s . t . x ∈ S \min f\left( x \right) \\ s.t. x\in S minf(x)s.t.x∈S 的局部最优解。
凸集与凸函数
定义设 f S ⊂ R n → R fS\subset R^{n}\rightarrow R fS⊂Rn→R, S S S是凸集如果对所有的 x ( 1 ) , x ( 2 ) ∈ S x^{(1)},x^{(2)}\in S x(1),x(2)∈S, λ ∈ ( 0 , 1 ) \lambda\in(0,1) λ∈(0,1),都有 f ( λ x ( 1 ) ( 1 − λ ) x ( 2 ) ) ⩽ λ f ( x ( 1 ) ) ( 1 − λ ) f ( x ( 2 ) ) f(\lambda x^{(1)}(1-\lambda )x^{(2)})\leqslant \lambda f( x^{(1)})(1-\lambda )f(x^{(2)}) f(λx(1)(1−λ)x(2))⩽λf(x(1))(1−λ)f(x(2)) 称 f ( x ) f(x) f(x)为 S S S上的凸函数。 Hessian矩阵 f f f存在二阶偏导数 x ∈ R n x\in R^{n} x∈Rn: ∇ 2 f ( x ) ( ∂ 2 f ( x ) ∂ x 1 2 ∂ 2 f ( x ) ∂ x 1 ∂ x 2 . . . ∂ 2 f ( x ) ∂ x 1 ∂ x n ∂ 2 f ( x ) ∂ x 2 ∂ x 1 ∂ 2 f ( x ) ∂ x 2 2 . . . ∂ 2 f ( x ) ∂ x 2 ∂ x n ⋮ ⋮ ⋱ ⋮ ∂ 2 f ( x ) ∂ x n ∂ x 1 ∂ 2 f ( x ) ∂ x n ∂ x 2 . . . ∂ 2 f ( x ) ∂ x n 2 ) \nabla ^2f\left( x \right) \left( \begin{matrix} \frac{\partial ^2f\left( x \right)}{\partial x_{1}^{2}} \frac{\partial ^2f\left( x \right)}{\partial x_1\partial x_2} ... \frac{\partial ^2f\left( x \right)}{\partial x_1\partial x_n}\\ \frac{\partial ^2f\left( x \right)}{\partial x_2\partial x_1} \frac{\partial ^2f\left( x \right)}{\partial x_{2}^{2}} ... \frac{\partial ^2f\left( x \right)}{\partial x_2\partial x_n}\\ \vdots \vdots \ddots \vdots\\ \frac{\partial ^2f\left( x \right)}{\partial x_n\partial x_1} \frac{\partial ^2f\left( x \right)}{\partial x_n\partial x_2} ... \frac{\partial ^2f\left( x \right)}{\partial x_{n}^{2}}\\ \end{matrix} \right) ∇2f(x) ∂x12∂2f(x)∂x2∂x1∂2f(x)⋮∂xn∂x1∂2f(x)∂x1∂x2∂2f(x)∂x22∂2f(x)⋮∂xn∂x2∂2f(x)......⋱...∂x1∂xn∂2f(x)∂x2∂xn∂2f(x)⋮∂xn2∂2f(x) 下面给出一系列定理
定理1有限凸函数的线性组合为凸函数定理2凸函数在有限函数值处的沿任意方向导数都存在定理3凸函数在定义区间上必连续定理4设 S S S是 R n R^{n} Rn上的非空凸集 f ( x ) f(x) f(x)是定义在 S S S上的凸函数则 f ( x ) f(x) f(x)在 S S S上的局限极小点是全局极小点且极小点的集合为凸集凸函数的一阶、二阶判断定理 ∃ x ( 1 ) ≠ x ( 2 ) ∈ S ≠ ⊘ , S 是一个凸集 f 可微 : f 为凸函数 ⇔ f ( x ( 2 ) ) ⩾ f ( x ( 1 ) ) ∇ f ( x ( 1 ) ) T ( x ( 2 ) − x ( 1 ) ) \exists x^{\left( 1 \right)}\ne x^{\left( 2 \right)}\in S\ne \oslash ,S\text{是一个凸集}f\text{可微}:f\text{为凸函数}\Leftrightarrow f\left( x^{\left( 2 \right)} \right) \geqslant f\left( x^{\left( 1 \right)} \right) \nabla f\left( x^{\left( 1 \right)} \right) ^T\left( x^{\left( 2 \right)}-x^{\left( 1 \right)} \right) ∃x(1)x(2)∈S⊘,S是一个凸集f可微:f为凸函数⇔f(x(2))⩾f(x(1))∇f(x(1))T(x(2)−x(1))在凸集定义内的点的Hessian矩阵半正定 ⇔ \Leftrightarrow ⇔凸函数在凸集定义上所有的点在其处的Hessian矩阵正定 ⇒ \Rightarrow ⇒严格凸函数 凸规划 min f ( x ) s . t . { g i ( x ) ⩾ 0 , i 1 , 2 , . . . , m h j ( x ) 0 , j 1 , 2 , . . . , n \min f\left( x \right) \\ s.t.\left\{ \begin{array}{c} g_i\left( x \right) \geqslant 0,i1,2,...,m\\ h_j\left( x \right) 0,j1,2,...,n\\ \end{array} \right. minf(x)s.t.{gi(x)⩾0,i1,2,...,mhj(x)0,j1,2,...,n f ( x ) f\left( x \right) f(x)为凸函数 g i ( x ) g_i\left( x \right) gi(x)为凹函数 h j ( x ) h_j\left( x \right) hj(x)为线性函数称此类问题为凸规划
极值条件
无约束条件的极值判断条件
无约束非线性规划的一般数学模型为如下形式 min f ( x ) x ∈ R n \min f\left( x \right) \\ x\in \mathbb{R} ^n minf(x)x∈Rn 定理1 设函数 f ( x ) 在 x ‾ 处可微 , 若存在向量 d , 使得 ∇ f ( x ‾ ) T d 0 , 则存在数 δ 0 , 使得对 ∀ λ ∈ ( 0 , δ ) , 有 f ( x ‾ λ d ) f ( x ‾ ) \text{设函数}f\left( x \right) \text{在}\overline{x}\text{处可微},\text{若存在向量}d,\text{使得}\nabla f\left( \overline{x} \right) ^Td0,\text{则存在数}\delta 0,\text{使得对}\forall \lambda \in \left( 0,\delta \right) ,\text{有}f\left( \overline{x}\lambda d \right) f\left( \overline{x} \right) 设函数f(x)在x处可微,若存在向量d,使得∇f(x)Td0,则存在数δ0,使得对∀λ∈(0,δ),有f(xλd)f(x) 定理2 设函数 f ( x ) 在点 x ‾ 处可微若 x ‾ 是局部最小点 , 则 ∇ f ( x ‾ ) 0 \text{设函数}f\left( x \right) \text{在点}\overline{x}\text{处可微若}\overline{x}\text{是局部最小点},\text{则}\nabla f\left( \overline{x} \right) 0 设函数f(x)在点x处可微若x是局部最小点,则∇f(x)0 定理3 设函数 f ( x ) 在点 x ‾ 处二次可微 , 若 x ‾ 是局部极小点则 ∇ f ( x ‾ ) 0 , 且 ∇ 2 f ( x ‾ ) 半正定 \text{设函数}f\left( x \right) \text{在点}\overline{x}\text{处二次可微},\text{若}\overline{x}\text{是局部极小点则}\nabla f\left( \overline{x} \right) 0,\text{且}\nabla ^2f\left( \overline{x} \right) \text{半正定} 设函数f(x)在点x处二次可微,若x是局部极小点则∇f(x)0,且∇2f(x)半正定 定理4——局部极小点 设函数 f ( x ) 在点 x ‾ 处二次可微 , 若 ∇ f ( x ‾ ) 0 , 且 ∇ 2 f ( x ‾ ) 正定则 x ‾ 是局部极小点 \text{设函数}f\left( x \right) \text{在点}\overline{x}\text{处二次可微},\text{若}\nabla f\left( \overline{x} \right) 0,\text{且}\nabla ^2f\left( \overline{x} \right) \text{正定则}\overline{x}\text{是局部极小点} 设函数f(x)在点x处二次可微,若∇f(x)0,且∇2f(x)正定则x是局部极小点
定理5——全局最小点 设函数 f ( x ) 在 R n 上的可微函数 , x ‾ ∈ R n , 则 x ‾ 是全局极小点 ⇔ ∇ f ( x ‾ ) 0 \text{设函数}f\left( x \right) \text{在}R^n\text{上的可微函数},\overline{x}\in R^n,\text{则}\overline{x}\text{是全局极小点}\Leftrightarrow \nabla f\left( \overline{x} \right) 0 设函数f(x)在Rn上的可微函数,x∈Rn,则x是全局极小点⇔∇f(x)0
有约束条件的极值判断条件
有约束非线性规划的一般数学模型为如下表现形式 min f ( x ) s . t . { g i ( x ) ⩾ 0 , i 1 , 2 , . . . , m h j ( x ) 0 , j 1 , 2 , . . . , n \min f\left( x \right) \\ s.t.\left\{ \begin{array}{c} g_i\left( x \right) \geqslant 0,i1,2,...,m\\ h_j\left( x \right) 0,j1,2,...,n\\ \end{array} \right. minf(x)s.t.{gi(x)⩾0,i1,2,...,mhj(x)0,j1,2,...,n 其中 x ∈ R n x\in R^n x∈Rn f ( x ) f(x) f(x) g i ( x ) g_i(x) gi(x) h j ( x ) h_j(x) hj(x)都是定义在 R n R^n Rn上的实值函数。 FRitz John条件——一般约束的一阶必要条件 x ‾ 为可行点 , I { i ∣ g i ( x ‾ ) 0 } , f 和 g i ( x ‾ ) ( i ∈ I ) 在点 x ‾ 处可微 , g i ( x ‾ ) ( i ∉ I ) 在点 x ‾ 处连续 h j ( x ) ( j 1 , 2 , . . . , l ) 在点 x ‾ 处连续可微。如果 x ‾ 是局部最优解 , 则存在不全为 0 的数 u 0 , u i ( i ∈ I ) 和 v j ( j 1 , 2 , . . . , l ) , 使得 u 0 ∇ f ( x ‾ ) ∑ i ∈ I u i ∇ g i ( x ∗ ) − ∑ j 1 l μ j ∇ h j ( x ∗ ) 0 , u 0 , u i ⩾ 0 \overline{x}\text{为可行点},I\left\{ i|g_i\left( \overline{x} \right) 0 \right\} ,f\text{和}g_i\left( \overline{x} \right) \left( i\in I \right) \text{在点}\overline{x}\text{处可微},g_i\left( \overline{x} \right) \left( i\notin I \right) \text{在点}\overline{x}\text{处连续}h_j\left( x \right) \left( j1,2,...,l \right) \text{在点}\overline{x}\text{处连续可微。如果}\overline{x}\text{是局部最优解},\text{则存在不全为}0\text{的数}u_0,u_i\left( i\in I \right) \text{和}v_j\left( j1,2,...,l \right) ,\text{使得} \\ u_0\nabla f(\overline{x})\sum_{i\in I}^{}{u_i\nabla g_i(x^*)}-\sum_{j1}^l{\mu _j\nabla h_j(x^*)}0, u_0,u_i\geqslant 0 x为可行点,I{i∣gi(x)0},f和gi(x)(i∈I)在点x处可微,gi(x)(i∈/I)在点x处连续hj(x)(j1,2,...,l)在点x处连续可微。如果x是局部最优解,则存在不全为0的数u0,ui(i∈I)和vj(j1,2,...,l),使得u0∇f(x)i∈I∑ui∇gi(x∗)−j1∑lμj∇hj(x∗)0,u0,ui⩾0 无约束非线性规划
一维搜索
一维搜索又称直线搜索是指单变量函数最优化即求一维问题的最优解方法这里介绍黄金分割法
先在区间 [ a , b ] \left[ a,b \right] [a,b]上插入两个试探点 t l a ( 1 − β ) ( b − a ) , t r a β ( b − a ) t_la\left( 1-\beta \right) \left( b-a \right) ,t_ra\beta \left( b-a \right) tla(1−β)(b−a),traβ(b−a)由单峰函数性质丢掉 1 3 \frac{1}{3} 31区间逐步搜索 黄金分割法是一维搜索中用于求解单变量函数最优解的一种高效方法以下是对它更详细的介绍
算法步骤
初始化 给定初始搜索区间 [ a , b ] [a,b] [a,b]以及精度要求 ϵ \epsilon ϵ。 计算试探点 计算两个试探点 t l a ( 1 − β ) ( b − a ) t_la(1 - \beta)(b - a) tla(1−β)(b−a) t r a β ( b − a ) t_ra\beta(b - a) traβ(b−a)。 比较函数值 计算 f ( t l ) f(t_l) f(tl)和 f ( t r ) f(t_r) f(tr)比较它们的大小。若 f ( t l ) f ( t r ) f(t_l)f(t_r) f(tl)f(tr)则新的搜索区间为 [ a , t r ] [a,t_r] [a,tr]。若 f ( t l ) ≥ f ( t r ) f(t_l)\geq f(t_r) f(tl)≥f(tr)则新的搜索区间为 [ t l , b ] [t_l,b] [tl,b]。 判断终止条件 检查新的搜索区间长度是否小于精度要求 ϵ \epsilon ϵ即 ∣ b − a ∣ ϵ |b - a|\epsilon ∣b−a∣ϵ。若是则取当前区间的中点或两个试探点中函数值较好的点作为近似最优解算法结束。若否则返回步骤2继续迭代。
示例
假设要求函数 f ( x ) − x 2 4 x − 3 f(x)-x^2 4x - 3 f(x)−x24x−3在区间 [ 0 , 4 ] [0,4] [0,4]上的最大值。
首先初始化 a 0 a 0 a0 b 4 b 4 b4 β 0.618 \beta 0.618 β0.618。计算试探点 t l 0 ( 1 − 0.618 ) × ( 4 − 0 ) 1.528 t_l0(1 - 0.618)\times(4 - 0)1.528 tl0(1−0.618)×(4−0)1.528。 t r 0 0.618 × ( 4 − 0 ) 2.472 t_r00.618\times(4 - 0)2.472 tr00.618×(4−0)2.472。 计算函数值 f ( t l ) f ( 1.528 ) − 1.52 8 2 4 × 1.528 − 3 1.315 f(t_l)f(1.528)-1.528^2 4\times1.528 - 31.315 f(tl)f(1.528)−1.52824×1.528−31.315。 f ( t r ) f ( 2.472 ) − 2.47 2 2 4 × 2.472 − 3 1.985 f(t_r)f(2.472)-2.472^2 4\times2.472 - 31.985 f(tr)f(2.472)−2.47224×2.472−31.985。 因为 f ( t l ) f ( t r ) f(t_l)f(t_r) f(tl)f(tr)所以新的搜索区间为 [ 0 , 2.472 ] [0,2.472] [0,2.472]。继续迭代直到满足精度要求。
特点
优点 不需要计算函数的导数适用于导数不存在或难以计算的函数。算法简单易于实现。收敛速度较快每次迭代都能将搜索区间缩小一定比例。 缺点 对于非单峰函数可能无法找到全局最优解。收敛速度相对较慢与一些基于导数的方法相比需要更多的迭代次数才能达到相同的精度。
代码模板
import math# 定义黄金分割比
GOLDEN_RATIO (math.sqrt(5) - 1) / 2def golden_section_search(func, a, b, epsilon1e-6):黄金分割法进行一维搜索:param func: 目标函数:param a: 搜索区间左端点:param b: 搜索区间右端点:param epsilon: 精度要求:return: 近似最优解# 计算初始试探点t_l a (1 - GOLDEN_RATIO) * (b - a)t_r a GOLDEN_RATIO * (b - a)# 计算初始试探点的函数值f_l func(t_l)f_r func(t_r)while (b - a) epsilon:if f_l f_r:# 新的搜索区间为 [a, t_r]b t_rt_r t_lf_r f_lt_l a (1 - GOLDEN_RATIO) * (b - a)f_l func(t_l)else:# 新的搜索区间为 [t_l, b]a t_lt_l t_rf_l f_rt_r a GOLDEN_RATIO * (b - a)f_r func(t_r)# 返回近似最优解return (a b) / 2# 示例目标函数
def example_function(x):return -x**2 4*x - 3# 测试代码
if __name__ __main__:a 0b 4result golden_section_search(example_function, a, b)print(f近似最优解: {result})print(f最优解处的函数值: {example_function(result)})最速下降法
最速下降法属于迭代算法。求解无约束非线性规划的迭代算法一般具有如下形式
选定初始点 x ( 0 ) , k 0 x^{(0)},k0 x(0),k0 k k 1 kk1 kk1寻找一个合适的方向 P ( k ) P^{(k)} P(k)求出沿着这个方向的前进步长 λ ( k ) \lambda^{(k)} λ(k)得到新的点 x ( k 1 ) x ( k ) λ ( k ) P ( k ) x^{(k1)}x^{(k)}\lambda^{(k)}P^{(k)} x(k1)x(k)λ(k)P(k) 最速下降算法的基本思想是沿着负梯度方向搜索极值点即搜索方向 P ( k ) − ∇ f ( x ( k ) ) P^{(k)}-\nabla f(x^{(k)}) P(k)−∇f(x(k)) 搜索步长满足以下条件: f ( x ( k ) λ ( k ) P ( k ) ) m i n λ ⩾ 0 f ( x ( k ) − λ ∇ f ( x ( k ) ) ) f(x^{(k)}\lambda ^{(k)}P^{(k)})\underset{\lambda\geqslant 0 }{min}f(x^{(k)}-\lambda \nabla f(x^{(k)})) f(x(k)λ(k)P(k))λ⩾0minf(x(k)−λ∇f(x(k)))
算法详细步骤
初始化 选定初始点 x ( 0 ) x^{(0)} x(0)并令迭代计数器 k 0 k 0 k0。初始点的选择通常会影响算法的收敛速度和最终结果在实际应用中可以根据问题的特点或经验来选择合适的初始点。迭代过程
步骤一更新迭代计数器 k k 1 kk 1 kk1 步骤二确定搜索方向 计算当前点 x ( k ) x^{(k)} x(k) 处目标函数 f ( x ) f(x) f(x) 的梯度 ∇ f ( x ( k ) ) \nabla f(x^{(k)}) ∇f(x(k))搜索方向 P ( k ) P^{(k)} P(k) 取为负梯度方向即 P ( k ) − ∇ f ( x ( k ) ) P^{(k)}-\nabla f(x^{(k)}) P(k)−∇f(x(k))。梯度 ∇ f ( x ) \nabla f(x) ∇f(x) 是一个向量它的每个分量是目标函数关于各个变量的偏导数。 步骤三确定搜索步长 搜索步长 λ ( k ) \lambda^{(k)} λ(k) 要满足 f ( x ( k ) λ ( k ) P ( k ) ) m i n λ ⩾ 0 f ( x ( k ) − λ ∇ f ( x ( k ) ) ) f(x^{(k)}\lambda ^{(k)}P^{(k)})\underset{\lambda\geqslant 0 }{min}f(x^{(k)}-\lambda \nabla f(x^{(k)})) f(x(k)λ(k)P(k))λ⩾0minf(x(k)−λ∇f(x(k)))。这意味着要在负梯度方向上找到一个最优的步长使得目标函数值下降最多。在实际计算中精确求解这个一维搜索问题可能比较困难因此常常采用一些近似的方法如黄金分割法、斐波那契法等进行一维搜索来确定步长。 步骤四更新迭代点 得到新的点 x ( k 1 ) x ( k ) λ ( k ) P ( k ) x^{(k 1)}x^{(k)}\lambda^{(k)}P^{(k)} x(k1)x(k)λ(k)P(k)。
终止条件 可以根据不同的准则来判断迭代是否终止常见的终止条件有
梯度的模小于某个给定的小正数 ϵ \epsilon ϵ即 ∥ ∇ f ( x ( k ) ) ∥ ϵ \left\|\nabla f(x^{(k)})\right\|\epsilon ∇f(x(k)) ϵ。这表示目标函数在当前点的变化率已经很小可能已经接近极值点。相邻两次迭代点的距离小于某个给定的小正数 ϵ \epsilon ϵ即 ∥ x ( k 1 ) − x ( k ) ∥ ϵ \left\|x^{(k 1)}-x^{(k)}\right\|\epsilon x(k1)−x(k) ϵ。目标函数值的变化小于某个给定的小正数 ϵ \epsilon ϵ即 ∣ f ( x ( k 1 ) ) − f ( x ( k ) ) ∣ ϵ \left|f(x^{(k 1)})-f(x^{(k)})\right|\epsilon f(x(k1))−f(x(k)) ϵ。
代码实现示例
import numpy as np
import math# 定义黄金分割比
GOLDEN_RATIO (math.sqrt(5) - 1) / 2# 黄金分割法进行一维搜索
def golden_section_search(func, a0, b10, epsilon1e-6):t_l a (1 - GOLDEN_RATIO) * (b - a)t_r a GOLDEN_RATIO * (b - a)f_l func(t_l)f_r func(t_r)while (b - a) epsilon:if f_l f_r:b t_rt_r t_lf_r f_lt_l a (1 - GOLDEN_RATIO) * (b - a)f_l func(t_l)else:a t_lt_l t_rf_l f_rt_r a GOLDEN_RATIO * (b - a)f_r func(t_r)return (a b) / 2# 最速下降法
def steepest_descent(f, grad_f, x0, epsilon1e-6, max_iter1000):最速下降法求解无约束非线性规划问题:param f: 目标函数:param grad_f: 目标函数的梯度函数:param x0: 初始点:param epsilon: 收敛精度:param max_iter: 最大迭代次数:return: 最优解x np.array(x0, dtypenp.float64)k 0while k max_iter:# 计算当前点的梯度grad grad_f(x)# 判断是否满足收敛条件if np.linalg.norm(grad) epsilon:break# 确定搜索方向为负梯度方向p -grad# 定义一维搜索的目标函数def one_dim_func(lambda_):return f(x lambda_ * p)# 使用黄金分割法确定步长lambda_ golden_section_search(one_dim_func)# 更新迭代点x x lambda_ * pk 1return x# 示例目标函数
def example_f(x):示例目标函数f(x) x1^2 x2^2return np.sum(np.square(x))# 示例目标函数的梯度
def example_grad_f(x):示例目标函数的梯度grad f(x) [2*x1, 2*x2]return 2 * x# 主程序
if __name__ __main__:# 初始点x0 [1.0, 1.0]# 调用最速下降法求解result steepest_descent(example_f, example_grad_f, x0)print(最优解:, result)print(最优值:, example_f(result))最优步长的求解
以下是使用黄金分割法和斐波那契法求步长 λ ( k ) \lambda^{(k)} λ(k)的具体方法
黄金分割法
原理 黄金分割法是基于黄金分割比例来逐步缩小搜索区间以逼近函数的最优解。对于一个单峰函数在搜索区间 [ a , b ] [a,b] [a,b]内通过在两个特定点 x 1 x_1 x1和 x 2 x_2 x2处计算函数值根据函数值的大小关系不断舍弃不包含最优解的区间从而使搜索区间不断缩小。这两个特定点 x 1 x_1 x1和 x 2 x_2 x2与区间端点的距离满足黄金分割比例即 x 1 − a b − a b − x 2 b − a 5 − 1 2 ≈ 0.618 \frac{x_1 - a}{b - a}\frac{b - x_2}{b - a}\frac{\sqrt{5}-1}{2}\approx0.618 b−ax1−ab−ab−x225 −1≈0.618。 计算步骤 首先确定初始搜索区间 [ a , b ] [a,b] [a,b]即确定步长 λ \lambda λ的取值范围保证目标函数 f ( x ( k ) λ P ( k ) ) f(x^{(k)}\lambda P^{(k)}) f(x(k)λP(k))在该区间内是单峰的。计算区间内的两个试探点 x 1 a 0.382 ( b − a ) x_1 a 0.382(b - a) x1a0.382(b−a) x 2 a 0.618 ( b − a ) x_2 a 0.618(b - a) x2a0.618(b−a)并计算对应的函数值 f 1 f ( x ( k ) x 1 P ( k ) ) f_1 f(x^{(k)}x_1 P^{(k)}) f1f(x(k)x1P(k)) f 2 f ( x ( k ) x 2 P ( k ) ) f_2 f(x^{(k)}x_2 P^{(k)}) f2f(x(k)x2P(k))。比较 f 1 f_1 f1和 f 2 f_2 f2的大小如果 f 1 f 2 f_1 f_2 f1f2则新的搜索区间为 [ a , x 2 ] [a,x_2] [a,x2]否则新的搜索区间为 [ x 1 , b ] [x_1,b] [x1,b]。重复上述步骤不断缩小搜索区间直到满足预设的精度要求此时的区间中点或其中一个端点即可作为近似的最优步长 λ ( k ) \lambda^{(k)} λ(k)。
斐波那契法
原理 斐波那契法也是一种用于一维搜索的区间收缩方法它与黄金分割法类似但在确定试探点的位置时是根据斐波那契数列来确定的。斐波那契数列满足 F n F n − 1 F n − 2 F_n F_{n - 1}F_{n - 2} FnFn−1Fn−2其中 F 0 0 F_0 0 F00 F 1 1 F_1 1 F11。在搜索过程中利用斐波那契数列的性质来确定试探点的位置使得在给定的搜索次数下能够最大程度地缩小搜索区间。 计算步骤 同样先确定初始搜索区间 [ a , b ] [a,b] [a,b]和搜索次数 n n n。根据斐波那契数列计算出第 n n n项 F n F_n Fn以及前一项 F n − 1 F_{n - 1} Fn−1第一个试探点为 x 1 a F n − 1 F n ( b − a ) x_1 a\frac{F_{n - 1}}{F_n}(b - a) x1aFnFn−1(b−a)第二个试探点为 x 2 a F n − 2 F n ( b − a ) x_2 a\frac{F_{n - 2}}{F_n}(b - a) x2aFnFn−2(b−a)计算对应的函数值 f 1 f ( x ( k ) x 1 P ( k ) ) f_1 f(x^{(k)}x_1 P^{(k)}) f1f(x(k)x1P(k)) f 2 f ( x ( k ) x 2 P ( k ) ) f_2 f(x^{(k)}x_2 P^{(k)}) f2f(x(k)x2P(k))。比较 f 1 f_1 f1和 f 2 f_2 f2的大小根据大小关系舍弃相应的区间更新搜索区间和试探点。随着搜索次数的增加不断重复上述步骤当达到预定的搜索次数 n n n时得到的搜索区间的中点或其中一个端点就是近似的最优步长 λ ( k ) \lambda^{(k)} λ(k)。
牛顿法
牛顿法基本思想以目标函数的一个二次函数作为其近似然后求出这个二次函数的极小值点从而该极小值点近似为原目标函数的一个局部极小值点。 设 f ( x ) f(x) f(x)是二次可微函数。假设上一步迭代到了 x ( k ) x^{(k)} x(k)。我们在点 x ( k ) x^{(k)} x(k)处 T a y l o r Taylor Taylor展开取二阶近似: f ( x ) ≈ Q ( x ) f ( x ( k ) ) ∇ f ( x ( k ) ) T ( x − x ( k ) ) 1 2 ( x − x ( k ) ) T ∇ 2 f ( x ( k ) ) ( x − x ( k ) ) f\left( x \right) \approx Q\left( x \right) f\left( x^{\left( k \right)} \right) \nabla f\left( x^{\left( k \right)} \right) ^T\left( x-x^{\left( k \right)} \right) \frac{1}{2}\left( x-x^{\left( k \right)} \right) ^T\nabla ^2f\left( x^{\left( k \right)} \right) \left( x-x^{\left( k \right)} \right) f(x)≈Q(x)f(x(k))∇f(x(k))T(x−x(k))21(x−x(k))T∇2f(x(k))(x−x(k)) g ( x ( k ) ) Δ ∇ f ( x ( k ) ) , H ( x ( k ) ) Δ ∇ 2 f ( x ( k ) ) g\left( x^{\left( k \right)} \right) \overset{\varDelta}{}\nabla f\left( x^{\left( k \right)} \right) ,H\left( x^{\left( k \right)} \right) \overset{\varDelta}{}\nabla ^2f\left( x^{\left( k \right)} \right) g(x(k))Δ∇f(x(k)),H(x(k))Δ∇2f(x(k)) 反解可得到 x x ( k ) − ∇ 2 f ( x ( k ) ) − 1 ∇ f ( x ( k ) ) xx^{\left( k \right)}-\nabla ^2f\left( x^{\left( k \right)} \right) ^{-1}\nabla f\left( x^{\left( k \right)} \right) xx(k)−∇2f(x(k))−1∇f(x(k)) 基于上面给出的信息我们得到牛顿算法 1给定初始点 x ( 1 ) ∈ R n x^{\left( 1 \right)}\in R^n x(1)∈Rn给定误差 ϵ 0 \epsilon0 ϵ0,令 k 1 k1 k1 2由 ∇ 2 f ( x ( k ) ) d ( k ) − ∇ f ( x ( k ) ) \nabla ^2f\left( x^{\left( k \right)} \right) d^{\left( k \right)}-\nabla f\left( x^{\left( k \right)} \right) ∇2f(x(k))d(k)−∇f(x(k))解线性方程组得到搜索方向 d ( k ) d^{(k)} d(k) 3置 x ( k 1 ) x ( k ) d ( k ) x^{\left( k1 \right)}x^{\left( k \right)}d^{\left( k \right)} x(k1)x(k)d(k)即 x ( k 1 ) x ( k ) − ∇ 2 f ( x ( k ) ) − 1 ∇ f ( x ( k ) ) x^{(k1)}x^{\left( k \right)}-\nabla ^2f\left( x^{\left( k \right)} \right) ^{-1}\nabla f\left( x^{\left( k \right)} \right) x(k1)x(k)−∇2f(x(k))−1∇f(x(k)) 4当 ∥ d ( k ) ∥ ϵ \left\| d^{\left( k \right)} \right\| \epsilon d(k) ϵ迭代终止否则 k k 1 kk1 kk1,跳转到 2 2 2
import numpy as npdef newton_method(f, grad_f, hess_f, x1, epsilon1e-6, max_iter100):牛顿算法实现参数:f: 目标函数grad_f: 目标函数的梯度hess_f: 目标函数的海森矩阵x1: 初始点numpy 数组epsilon: 误差阈值默认为 1e-6max_iter: 最大迭代次数默认为 100返回:x: 最优解iter_count: 迭代次数x x1.copy()k 0while k max_iter:# 步骤 2: 求解线性方程组得到搜索方向grad grad_f(x)hess hess_f(x)try:# 解线性方程组 H * d -gradd np.linalg.solve(hess, -grad)except np.linalg.LinAlgError:print(海森矩阵不可逆无法求解线性方程组。)break# 步骤 3: 更新 xx_new x d# 步骤 4: 判断是否满足终止条件if np.linalg.norm(d) epsilon:breakx x_newk 1return x, k# 示例定义一个简单的目标函数及其梯度和海森矩阵
def f(x):# 目标函数 f(x) x1^2 x2^2return x[0]**2 x[1]**2def grad_f(x):# 目标函数的梯度return np.array([2 * x[0], 2 * x[1]])def hess_f(x):# 目标函数的海森矩阵return np.array([[2, 0], [0, 2]])# 初始点
x1 np.array([1.0, 1.0])# 调用牛顿算法
optimal_x, iter_count newton_method(f, grad_f, hess_f, x1)print(最优解:, optimal_x)
print(迭代次数:, iter_count)阻尼牛顿法
阻尼牛顿法与原始牛顿法的主要区别在于增加了沿着牛顿方向进行一维搜索其迭代公式: x ( k 1 ) x ( k ) λ ( k ) d ( k ) x^{\left( k1 \right)}x^{\left( k \right)}\lambda ^{\left( k \right)}d^{\left( k \right)} x(k1)x(k)λ(k)d(k) 1给定初始点 x ( 1 ) ∈ R n x^{\left( 1 \right)}\in R^n x(1)∈Rn给定误差 ϵ 0 \epsilon0 ϵ0,令 k 1 k1 k1 2由 ∇ 2 f ( x ( k ) ) d ( k ) − ∇ f ( x ( k ) ) \nabla ^2f\left( x^{\left( k \right)} \right) d^{\left( k \right)}-\nabla f\left( x^{\left( k \right)} \right) ∇2f(x(k))d(k)−∇f(x(k))解线性方程组得到搜索方向 d ( k ) d^{(k)} d(k) 3 从 x ( k ) x^{(k)} x(k)出发沿着 d ( k ) d^{(k)} d(k)方向一维搜索 min λ f ( x ( k ) λ d ( k ) ) \underset{\lambda}{\min}\,\,f\left( x^{\left( k \right)}\lambda d^{\left( k \right)} \right) λminf(x(k)λd(k)) 4置 x ( k 1 ) x ( k ) d ( k ) x^{\left( k1 \right)}x^{\left( k \right)}d^{\left( k \right)} x(k1)x(k)d(k)即 x ( k 1 ) x ( k ) − λ ( k ) ∇ 2 f ( x ( k ) ) − 1 ∇ f ( x ( k ) ) x^{(k1)}x^{\left( k \right)}-\lambda ^{\left( k \right)}\nabla ^2f\left( x^{\left( k \right)} \right) ^{-1}\nabla f\left( x^{\left( k \right)} \right) x(k1)x(k)−λ(k)∇2f(x(k))−1∇f(x(k)) 5当 ∥ d ( k ) ∥ ϵ \left\| d^{\left( k \right)} \right\| \epsilon d(k) ϵ迭代终止否则 k k 1 kk1 kk1,跳转到 2 2 2
模式搜索法
Powell方法
Powell方法是研究对称正定的二次函数 f ( x ) 1 2 x T Q x b T x c f(x)\frac{1}{2}x^TQxb^Txc f(x)21xTQxbTxc 极小化问题时候形成的。
1选定初始点 x ( 0 ) x^{\left( 0 \right)} x(0) n n n个线性无关向量 d ( i ) d^{(i)} d(i),给定误差 ϵ 0 \epsilon0 ϵ0,令 k 1 k1 k1, i 1 , 2 , . . . , n i1,2,...,n i1,2,...,n2依次对目标函数作直线搜索 { f ( x ( i − 1 ) λ i d ( i ) ) min λ ∈ R f ( x ( i − 1 ) λ d ( i ) ) x ( i ) x ( i − 1 ) λ i d ( i ) \left\{ \begin{array}{c} f\left( x^{\left( i-1 \right)}\lambda _id^{\left( i \right)} \right) \underset{\lambda \in R}{\min}f\left( x^{\left( i-1 \right)}\lambda d^{\left( i \right)} \right)\\ x^{\left( i \right)}x^{\left( i-1 \right)}\lambda _id^{\left( i \right)}\\ \end{array} \right. {f(x(i−1)λid(i))λ∈Rminf(x(i−1)λd(i))x(i)x(i−1)λid(i)3 对 i 1 , 2 , . . . , n − 1 i1,2,...,n-1 i1,2,...,n−1 d ( i ) d ( i 1 ) d^{(i)}d^{(i1)} d(i)d(i1), d ( n 1 ) x ( n ) − x ( 0 ) d^{(n1)}x^{(n)}-x^{(0)} d(n1)x(n)−x(0),置 d ( n ) d ( n 1 ) d^{(n)}d^{(n1)} d(n)d(n1)4直线搜索 x ( n 1 ) x ( n ) λ n 1 d ( n 1 ) x^{\left( n1 \right)}x^{\left( n \right)}\lambda _{n1}d^{\left( n1 \right)} x(n1)x(n)λn1d(n1) f ( x ( n ) λ n 1 d ( n 1 ) ) min f ( x ( n ) λ d ( n 1 ) ) f\left( x^{\left( n \right)}\lambda _{n1}d^{\left( n1 \right)} \right) \min f\left( x^{\left( n \right)}\lambda d^{\left( n1 \right)} \right) f(x(n)λn1d(n1))minf(x(n)λd(n1))5 ∥ x ( n 1 ) − x ( 0 ) ∥ ϵ , x ( n 1 ) \left\| x^{\left( n1 \right)}-x^{\left( 0 \right)} \right\| \epsilon ,x^{\left( n1 \right)} x(n1)−x(0) ϵ,x(n1)就是所求的极小值点则停止计算否则转第 6 6 6步6置 x ( 0 ) x ( n 1 ) x^{(0)}x^{(n1)} x(0)x(n1)转 ( 2 ) (2) (2)
有约束非线性规划
有约束非线性规划的一般数学模型为如下形式 min f ( x ) s . t . { g i ( x ) ⩾ 0 , i 1 , 2 , . . . , m h j ( x ) 0 , j 1 , 2 , . . . , n \min f\left( x \right) \\ s.t.\left\{ \begin{array}{c} g_i\left( x \right) \geqslant 0,i1,2,...,m\\ h_j\left( x \right) 0,j1,2,...,n\\ \end{array} \right. minf(x)s.t.{gi(x)⩾0,i1,2,...,mhj(x)0,j1,2,...,n 其中 f ( x ) g i ( x ) , h i ( x ) f(x)g_{i}(x),h_{i}(x) f(x)gi(x),hi(x)至少有一个非线性函数。
外点罚函数法
外点罚函数法是一种将有约束优化问题转化为无约束优化问题的方法。其核心思想是通过在目标函数中添加一个罚项当约束条件不满足时罚项的值会变得很大从而迫使无约束优化问题的解逐渐逼近原约束优化问题的可行域。
对于只含等式约束的优化问题 min f ( x ) s . t . h j ( x ) 0 , j 1 , 2 , . . . , l \min f\left( x \right) \\ s.t. h_j\left( x \right) 0,j1,2,...,l\\ minf(x)s.t.hj(x)0,j1,2,...,l 采用如下方法建立罚函数 F ( x , m k ) f ( x ) m k ∑ j 1 l h j 2 ( x ) F\left( x,m_k \right) f\left( x \right) m_k\sum_{j1}^l{h_{j}^{2}\left( x \right)} F(x,mk)f(x)mkj1∑lhj2(x) 这样我们就把原约束问题化为如下的无约束问题 m i n F ( x , m k ) min F(x,m_k) minF(x,mk)对于不等式条件 构造罚函数 F ( x , m k ) f ( x ) m k ∑ i 1 m g i 2 ( x ) u i ( g i ) , u 为单位跃阶函数 F\left( x,m_k \right) f\left( x \right) m_k\sum_{i1}^m{g_{i}^{2}\left( x \right) u_i\left( g_i \right)},u\text{为单位跃阶函数} F(x,mk)f(x)mki1∑mgi2(x)ui(gi),u为单位跃阶函数综合以上思想我们可以处理一般的罚函数 对于包含等式约束 h j ( x ) 0 h_j(x) 0 hj(x)0 j 1 , 2 , ⋯ , l j 1, 2, \cdots, l j1,2,⋯,l和不等式约束 g i ( x ) ≤ 0 g_i(x) \leq 0 gi(x)≤0 i 1 , 2 , ⋯ , m i 1, 2, \cdots, m i1,2,⋯,m的一般约束优化问题 min f ( x ) s.t. h j ( x ) 0 , j 1 , 2 , ⋯ , l g i ( x ) ≤ 0 , i 1 , 2 , ⋯ , m \begin{align*} \min \quad f(x) \\ \text{s.t.} \quad h_j(x) 0, \quad j 1, 2, \cdots, l \\ \quad g_i(x) \leq 0, \quad i 1, 2, \cdots, m \end{align*} mins.t.f(x)hj(x)0,j1,2,⋯,lgi(x)≤0,i1,2,⋯,m 可以构造外点罚函数 F ( x , m k ) f ( x ) m k ( ∑ j 1 l h j 2 ( x ) ∑ i 1 m [ max { 0 , g i ( x ) } ] 2 ) F(x, m_k) f(x) m_k \left( \sum_{j 1}^l h_j^2(x) \sum_{i 1}^m [\max\{0, g_i(x)\}]^2 \right) F(x,mk)f(x)mk(j1∑lhj2(x)i1∑m[max{0,gi(x)}]2) 其中 m k m_k mk 是罚因子且满足 0 m 1 m 2 ⋯ m k ⋯ 0 m_1 m_2 \cdots m_k \cdots 0m1m2⋯mk⋯ lim k → ∞ m k ∞ \lim_{k \to \infty} m_k \infty limk→∞mk∞。
算法步骤
初始化选择初始点 x ( 0 ) x^{(0)} x(0)初始罚因子 m 1 0 m_1 0 m10罚因子增长系数 c 1 c 1 c1收敛精度 ϵ 0 \epsilon 0 ϵ0令 k 1 k 1 k1。求解无约束优化问题以 x ( k − 1 ) x^{(k - 1)} x(k−1) 为初始点求解无约束优化问题 min F ( x , m k ) \min F(x, m_k) minF(x,mk)得到最优解 x ( k ) x^{(k)} x(k)。判断收敛条件如果满足 ∣ F ( x ( k ) , m k ) − F ( x ( k − 1 ) , m k − 1 ) ∣ ϵ \left| F(x^{(k)}, m_k) - F(x^{(k - 1)}, m_{k - 1}) \right| \epsilon F(x(k),mk)−F(x(k−1),mk−1) ϵ 或者约束条件的违反程度足够小例如 ∑ j 1 l h j 2 ( x ( k ) ) ∑ i 1 m [ max { 0 , g i ( x ( k ) ) } ] 2 ϵ \sum_{j 1}^l h_j^2(x^{(k)}) \sum_{i 1}^m [\max\{0, g_i(x^{(k)})\}]^2 \epsilon ∑j1lhj2(x(k))∑i1m[max{0,gi(x(k))}]2ϵ则停止迭代 x ( k ) x^{(k)} x(k) 即为原约束优化问题的近似最优解否则进入下一步。更新罚因子令 m k 1 c m k m_{k 1} c m_k mk1cmk k k 1 k k 1 kk1返回步骤 2。
Python 代码实现
import numpy as np
from scipy.optimize import minimizedef penalty_function(x, f, h, g, mk):计算外点罚函数的值:param x: 当前点:param f: 目标函数:param h: 等式约束函数列表:param g: 不等式约束函数列表:param mk: 罚因子:return: 罚函数的值penalty_term 0# 等式约束的罚项for hj in h:penalty_term hj(x) ** 2# 不等式约束的罚项for gi in g:penalty_term max(0, gi(x)) ** 2return f(x) mk * penalty_termdef outer_penalty_method(f, h, g, x0, m11.0, c10.0, epsilon1e-6, max_iter100):外点罚函数法求解约束优化问题:param f: 目标函数:param h: 等式约束函数列表:param g: 不等式约束函数列表:param x0: 初始点:param m1: 初始罚因子:param c: 罚因子增长系数:param epsilon: 收敛精度:param max_iter: 最大迭代次数:return: 最优解mk m1x np.array(x0)for k in range(max_iter):# 定义当前罚函数def current_penalty(x):return penalty_function(x, f, h, g, mk)# 求解无约束优化问题result minimize(current_penalty, x)x_new result.x# 判断收敛条件if k 0 and np.abs(current_penalty(x_new) - current_penalty(x)) epsilon:breakx x_new# 更新罚因子mk * creturn x# 示例求解一个简单的约束优化问题
# 目标函数
def f(x):return (x[0] - 1) ** 2 (x[1] - 2.5) ** 2# 等式约束
def h1(x):return x[0] x[1] - 3# 不等式约束
def g1(x):return -x[0] 2def g2(x):return -x[1] 2# 初始点
x0 [0, 0]# 调用外点罚函数法
optimal_x outer_penalty_method(f, [h1], [g1, g2], x0)print(最优解:, optimal_x)注意事项
罚因子的增长系数 c c c 和初始罚因子 m 1 m_1 m1 的选择会影响算法的收敛速度和稳定性需要根据具体问题进行调整。无约束优化问题的求解可以使用不同的优化算法这里使用了 scipy.optimize.minimize 函数你可以根据需要选择其他合适的算法。
内点罚函数法
内点罚函数法Interior Point Penalty Method也称为障碍函数法是一种用于求解有约束优化问题的重要方法。它主要适用于只包含不等式约束的优化问题基本思想是在可行域内部设置障碍使得迭代点始终保持在可行域内并且随着迭代的进行障碍的作用逐渐减弱从而使迭代点逼近最优解。以下将详细介绍内点罚函数法的原理、算法步骤并给出 Python 代码示例。
原理
考虑如下只含不等式约束的优化问题 min f ( x ) s.t. g i ( x ) ≤ 0 , i 1 , 2 , ⋯ , m \begin{align*} \min \quad f(x) \\ \text{s.t.} \quad g_i(x) \leq 0, \quad i 1, 2, \cdots, m \end{align*} mins.t.f(x)gi(x)≤0,i1,2,⋯,m 内点罚函数法通过构造障碍函数来实现。常见的障碍函数有倒数障碍函数和对数障碍函数。
倒数障碍函数 B ( x , r k ) f ( x ) r k ∑ i 1 m 1 − g i ( x ) B(x, r_k) f(x) r_k \sum_{i 1}^{m} \frac{1}{-g_i(x)} B(x,rk)f(x)rki1∑m−gi(x)1
对数障碍函数 B ( x , r k ) f ( x ) − r k ∑ i 1 m ln ( − g i ( x ) ) B(x, r_k) f(x) - r_k \sum_{i 1}^{m} \ln(-g_i(x)) B(x,rk)f(x)−rki1∑mln(−gi(x))
其中 r k r_k rk 是障碍因子且满足 r 1 r 2 ⋯ r k ⋯ r_1 r_2 \cdots r_k \cdots r1r2⋯rk⋯ lim k → ∞ r k 0 \lim_{k \to \infty} r_k 0 limk→∞rk0。当迭代点接近可行域边界时障碍函数的值会变得很大从而阻止迭代点越过边界。
算法步骤
初始化选择一个初始可行点 x ( 0 ) x^{(0)} x(0)即满足所有不等式约束 g i ( x ( 0 ) ) ≤ 0 g_i(x^{(0)}) \leq 0 gi(x(0))≤0 i 1 , 2 , ⋯ , m i 1, 2, \cdots, m i1,2,⋯,m初始障碍因子 r 1 0 r_1 0 r10障碍因子缩减系数 c ∈ ( 0 , 1 ) c \in (0, 1) c∈(0,1)收敛精度 ϵ 0 \epsilon 0 ϵ0令 k 1 k 1 k1。求解无约束优化问题以 x ( k − 1 ) x^{(k - 1)} x(k−1) 为初始点求解无约束优化问题 min B ( x , r k ) \min B(x, r_k) minB(x,rk)得到最优解 x ( k ) x^{(k)} x(k)。判断收敛条件如果满足 ∣ r k ∑ i 1 m 1 − g i ( x ( k ) ) ∣ ϵ \left| r_k \sum_{i 1}^{m} \frac{1}{-g_i(x^{(k)})} \right| \epsilon rk∑i1m−gi(x(k))1 ϵ对于倒数障碍函数或者 ∣ − r k ∑ i 1 m ln ( − g i ( x ( k ) ) ∣ ϵ \left| -r_k \sum_{i 1}^{m} \ln(-g_i(x^{(k)}) \right| \epsilon −rk∑i1mln(−gi(x(k)) ϵ对于对数障碍函数则停止迭代 x ( k ) x^{(k)} x(k) 即为原约束优化问题的近似最优解否则进入下一步。更新障碍因子令 r k 1 c r k r_{k 1} c r_k rk1crk k k 1 k k 1 kk1返回步骤 2。
Python 代码示例使用对数障碍函数
import numpy as np
from scipy.optimize import minimizedef logarithmic_barrier_function(x, f, g, rk):计算对数障碍函数的值:param x: 当前点:param f: 目标函数:param g: 不等式约束函数列表:param rk: 障碍因子:return: 对数障碍函数的值barrier_term 0for gi in g:if gi(x) 0:# 如果不满足不等式约束返回一个很大的值return np.infbarrier_term np.log(-gi(x))return f(x) - rk * barrier_termdef interior_point_method(f, g, x0, r11.0, c0.1, epsilon1e-6, max_iter100):内点罚函数法求解约束优化问题:param f: 目标函数:param g: 不等式约束函数列表:param x0: 初始可行点:param r1: 初始障碍因子:param c: 障碍因子缩减系数:param epsilon: 收敛精度:param max_iter: 最大迭代次数:return: 最优解rk r1x np.array(x0)for k in range(max_iter):# 定义当前对数障碍函数def current_barrier(x):return logarithmic_barrier_function(x, f, g, rk)# 求解无约束优化问题result minimize(current_barrier, x)x_new result.x# 判断收敛条件barrier_term 0for gi in g:barrier_term np.log(-gi(x_new))if np.abs(-rk * barrier_term) epsilon:breakx x_new# 更新障碍因子rk * creturn x# 示例求解一个简单的约束优化问题
# 目标函数
def f(x):return (x[0] - 1) ** 2 (x[1] - 2.5) ** 2# 不等式约束
def g1(x):return -x[0] 2def g2(x):return -x[1] 2# 初始可行点
x0 [1, 1]# 调用内点罚函数法
optimal_x interior_point_method(f, [g1, g2], x0)print(最优解:, optimal_x)代码说明
logarithmic_barrier_function 函数用于计算对数障碍函数的值。如果迭代点不满足不等式约束返回一个很大的值以避免该点被选为最优解。interior_point_method 函数实现了内点罚函数法的主要步骤包括初始化、求解无约束优化问题、判断收敛条件和更新障碍因子。示例定义了一个简单的约束优化问题包含两个不等式约束调用 interior_point_method 函数求解最优解。
注意事项
初始可行点的选择非常重要如果选择不当可能会导致算法无法收敛。障碍因子的缩减系数 c c c 和初始障碍因子 r 1 r_1 r1 的选择会影响算法的收敛速度和稳定性需要根据具体问题进行调整。 文章转载自： http://www.morning.kzcfp.cn.gov.cn.kzcfp.cn http://www.morning.nhgkm.cn.gov.cn.nhgkm.cn http://www.morning.mfmrg.cn.gov.cn.mfmrg.cn http://www.morning.qdlnw.cn.gov.cn.qdlnw.cn http://www.morning.srnth.cn.gov.cn.srnth.cn http://www.morning.ghslr.cn.gov.cn.ghslr.cn http://www.morning.yesidu.com.gov.cn.yesidu.com http://www.morning.mbfkt.cn.gov.cn.mbfkt.cn http://www.morning.kehejia.com.gov.cn.kehejia.com http://www.morning.rsnd.cn.gov.cn.rsnd.cn http://www.morning.chfxz.cn.gov.cn.chfxz.cn http://www.morning.pmlgr.cn.gov.cn.pmlgr.cn http://www.morning.ypzr.cn.gov.cn.ypzr.cn http://www.morning.tynqy.cn.gov.cn.tynqy.cn http://www.morning.nxtgb.cn.gov.cn.nxtgb.cn http://www.morning.dnhdp.cn.gov.cn.dnhdp.cn http://www.morning.drnjn.cn.gov.cn.drnjn.cn http://www.morning.dlwzm.cn.gov.cn.dlwzm.cn http://www.morning.nqcts.cn.gov.cn.nqcts.cn http://www.morning.wplbs.cn.gov.cn.wplbs.cn http://www.morning.clwhf.cn.gov.cn.clwhf.cn http://www.morning.tkxyx.cn.gov.cn.tkxyx.cn http://www.morning.ypklb.cn.gov.cn.ypklb.cn http://www.morning.wchcx.cn.gov.cn.wchcx.cn http://www.morning.rui931.cn.gov.cn.rui931.cn http://www.morning.xhklb.cn.gov.cn.xhklb.cn http://www.morning.nstml.cn.gov.cn.nstml.cn http://www.morning.xclgf.cn.gov.cn.xclgf.cn http://www.morning.1000sh.com.gov.cn.1000sh.com http://www.morning.wfkbk.cn.gov.cn.wfkbk.cn http://www.morning.xrpjr.cn.gov.cn.xrpjr.cn http://www.morning.swimstaracademy.cn.gov.cn.swimstaracademy.cn http://www.morning.ljcjc.cn.gov.cn.ljcjc.cn http://www.morning.dbfwq.cn.gov.cn.dbfwq.cn http://www.morning.xpzgg.cn.gov.cn.xpzgg.cn http://www.morning.bswnf.cn.gov.cn.bswnf.cn http://www.morning.bxfy.cn.gov.cn.bxfy.cn http://www.morning.lkrmp.cn.gov.cn.lkrmp.cn http://www.morning.hhpkb.cn.gov.cn.hhpkb.cn http://www.morning.ngqdp.cn.gov.cn.ngqdp.cn http://www.morning.rtbx.cn.gov.cn.rtbx.cn http://www.morning.jppb.cn.gov.cn.jppb.cn http://www.morning.fpxsd.cn.gov.cn.fpxsd.cn http://www.morning.xppj.cn.gov.cn.xppj.cn http://www.morning.rbbgh.cn.gov.cn.rbbgh.cn http://www.morning.jkcpl.cn.gov.cn.jkcpl.cn http://www.morning.pdbgm.cn.gov.cn.pdbgm.cn http://www.morning.wkrkb.cn.gov.cn.wkrkb.cn http://www.morning.nrjr.cn.gov.cn.nrjr.cn http://www.morning.wrwcf.cn.gov.cn.wrwcf.cn http://www.morning.gqfbh.cn.gov.cn.gqfbh.cn http://www.morning.c-ae.cn.gov.cn.c-ae.cn http://www.morning.wflsk.cn.gov.cn.wflsk.cn http://www.morning.rkck.cn.gov.cn.rkck.cn http://www.morning.bpmnh.cn.gov.cn.bpmnh.cn http://www.morning.qctsd.cn.gov.cn.qctsd.cn http://www.morning.rshkh.cn.gov.cn.rshkh.cn http://www.morning.nynpf.cn.gov.cn.nynpf.cn http://www.morning.ljdtn.cn.gov.cn.ljdtn.cn http://www.morning.jkzjs.cn.gov.cn.jkzjs.cn http://www.morning.mzkn.cn.gov.cn.mzkn.cn http://www.morning.srckl.cn.gov.cn.srckl.cn http://www.morning.jqswf.cn.gov.cn.jqswf.cn http://www.morning.kpyyf.cn.gov.cn.kpyyf.cn http://www.morning.thpzn.cn.gov.cn.thpzn.cn http://www.morning.gfkb.cn.gov.cn.gfkb.cn http://www.morning.mhnb.cn.gov.cn.mhnb.cn http://www.morning.wrkhf.cn.gov.cn.wrkhf.cn http://www.morning.cctgww.cn.gov.cn.cctgww.cn http://www.morning.hfyll.cn.gov.cn.hfyll.cn http://www.morning.mljtx.cn.gov.cn.mljtx.cn http://www.morning.trsdm.cn.gov.cn.trsdm.cn http://www.morning.mswkd.cn.gov.cn.mswkd.cn http://www.morning.gywfp.cn.gov.cn.gywfp.cn http://www.morning.pcrzf.cn.gov.cn.pcrzf.cn http://www.morning.cfnht.cn.gov.cn.cfnht.cn http://www.morning.lpzyq.cn.gov.cn.lpzyq.cn http://www.morning.txtgy.cn.gov.cn.txtgy.cn http://www.morning.yjprj.cn.gov.cn.yjprj.cn http://www.morning.kzxlc.cn.gov.cn.kzxlc.cn
