【最优化笔记3】线性规划 您所在的位置:网站首页 对偶单纯形法matlab程序 【最优化笔记3】线性规划

【最优化笔记3】线性规划

2024-07-01 16:47| 来源: 网络整理| 查看: 265

单纯形法是求解线性规划问题的一种通用的有效算法(必考点)。 此篇博客以我的好友Cizeron总结为基础完成,特此表示感谢。 在这里插入图片描述

目录 1.前置概念2.基本思想3.算法步骤4.算例5.算法收敛性6.Matlab实现1.输入问题2.建立初始单纯形表3.迭代寻找最优解4.输出结果5.附录(代码总表)

1.前置概念

(1)约束方程的规范形式: { m i n   f = c T x s . t .   A x = b , x ≥ 0   (1) \begin{cases}min \,f=c^Tx \\ s.t.\,Ax=b,& \text{$x\geq0 $ } \end{cases} \tag{1} {minf=cTxs.t.Ax=b,​x≥0 ​(1) 线性规划(1)的约束条件系数矩阵A通过初等行变换,总可以化为 [ I m , N ] [I_m,N] [Im​,N] 则约束条件可以写为 [ I m , N ] [ x B x N ] = b ′ (2) [I_m,N] \begin{bmatrix} x_B\\ x_N\\ \end{bmatrix}=b' \tag{2} [Im​,N][xB​xN​​]=b′(2) 其中,显然 B = I m B=I_m B=Im​即为基, x = ( x B = b ′ , 0 ) T x=(x_B=b',0)^T x=(xB​=b′,0)T为(1)的关于B基本解。将线性规划问题化为此规范形式是下文编程实现的第一步。 (2)基变换: 若初始基本可行解 x ( 0 ) x^{(0)} x(0)不是最优解,那么就还要找一个新的基本可行解:从原可行基中换出一个列向量(离基变量),再换入一个新的列向量(进基变量)(要保证线性无关),从而得到一个新的可行基,这就是基变换。

2.基本思想

单纯形法的基本思想是:给出一种规则,使由 LP问题一个基本可行解(极点)转移到另一个基本可行解,目标函数值是减小的,而且两个基本可行解之间的转换是容易实现的,经过有限次迭代,即可求得所需的最优基本可行解。

再具体一点来说: 对于一个优化问题 { m i n   f = c T x s . t .   A x = b , x ≥ 0   (1) \begin{cases}min \,f=c^Tx \\ s.t.\,Ax=b,& \text{$x\geq0 $ } \end{cases} \tag{1} {minf=cTxs.t.Ax=b,​x≥0 ​(1) 这里 b ≥ 0 , r a n k ( A ) = m b\geq 0,rank(A)=m b≥0,rank(A)=m,记 A = [ B , N ] = [ a 11 a 12 . . . a 1 n a 21 a 22 . . . a 2 n . . . . . . . . . . . . a m 1 a m 2 . . . a m n ] = ( p 1 , p 2 , . . . , p n ) A=[B,N]=\begin{bmatrix} a_{11} & a_{12}&...& a_{1n} \\ a_{21} & a_{22}&...& a_{2n} \\ ... &...&...& ...\\a_{m1} & a_{m2}&...& a_{mn}\end{bmatrix}=(p_1,p_2,...,p_n) A=[B,N]=⎣⎢⎢⎡​a11​a21​...am1​​a12​a22​...am2​​............​a1n​a2n​...amn​​⎦⎥⎥⎤​=(p1​,p2​,...,pn​),其中 r a n k ( B ) = m rank(B)=m rank(B)=m,即B为初始可行基。 设 x ( 0 ) = [ B − 1 b 0 ] x^{(0)}=\begin{bmatrix} B^{-1}b \\ 0 \end{bmatrix} x(0)=[B−1b0​]为问题(1)的一个初始可行解,则在 x ( 0 ) x^{(0)} x(0)处的目标函数值 f 0 = c T x ( 0 ) = [ c B , c N ] [ B − 1 b 0 ] = c B B − 1 b f_0=c^Tx^{(0)}=[c_B,c_N]\begin{bmatrix} B^{-1}b \\ 0 \end{bmatrix}=c_BB^{-1}b f0​=cTx(0)=[cB​,cN​][B−1b0​]=cB​B−1b。 若初始可行解 x ( 0 ) x^{(0)} x(0)不是最优解,则需要找到一个新的基本可行解,即进行基变换。 设 x = [ x B x N ] x=\begin{bmatrix} x_B \\ x_N \end{bmatrix} x=[xB​xN​​]为任意一个可行解,则由 A x = b Ax=b Ax=b得: x B = B − 1 b − B − 1 N x N x_B=B^{-1}b-B^{-1}Nx_N xB​=B−1b−B−1NxN​,在该 x x x点的目标函数值为: f = c T x = c B x B + c N x N = c B ( B − 1 b − B − 1 N x N ) + c N x N = c B B − 1 b − ( c B B − 1 N − c N ) x N = c B ( B − 1 b − B − 1 N x N ) + c N x N = c B B − 1 b − ∑ j ∈ R N ( c B B − 1 p j − c j ) x j = f 0 − ∑ j ∈ R N ( z j − c j ) x j f=c^Tx=c_Bx_B+c_Nx_N=c_B(B^{-1}b-B^{-1}Nx_N)+c_Nx_N=c_BB^{-1}b-(c_BB^{-1}N-c_N)x_N=c_B(B^{-1}b-B^{-1}Nx_N)+c_Nx_N=c_BB^{-1}b-\sum_{j\in R_N}(c_BB^{-1}p_j-c_j)x_j=f_0-\sum_{j\in R_N}(z_j-c_j)x_j f=cTx=cB​xB​+cN​xN​=cB​(B−1b−B−1NxN​)+cN​xN​=cB​B−1b−(cB​B−1N−cN​)xN​=cB​(B−1b−B−1NxN​)+cN​xN​=cB​B−1b−∑j∈RN​​(cB​B−1pj​−cj​)xj​=f0​−∑j∈RN​​(zj​−cj​)xj​,其中 R N R_N RN​为非基变量的下标集以及 z = c B B − 1 p j z=c_BB^{-1}p_j z=cB​B−1pj​。最后的可行解函数值为: f = f 0 − ∑ j ∈ R N ( z j − c j ) x j f=f_0-\sum_{j\in R_N}(z_j-c_j)x_j f=f0​−j∈RN​∑​(zj​−cj​)xj​因而,适当选取非基变量 x j ( j ∈ R N ) x_j(j\in R_N) xj​(j∈RN​)的数值就有可能使得 ∑ j ∈ R N ( z j − c j ) x j > 0 \sum_{j\in R_N}(z_j-c_j)x_j >0 ∑j∈RN​​(zj​−cj​)xj​>0从而达到目标函数下降的目的。为此,我们在原来的 n − m n-m n−m个非基变量中,令 n − m − 1 n-m-1 n−m−1个变量仍取0值,而那个“天选之子”(比如 x k x_k xk​)取正值,只要 ( z k − c k ) > 0 (z_k-c_k)>0 (zk​−ck​)>0就能达到优化的目的。 至于这个 k k k怎么选,怎么生成新的基本可行解,可以看看上一节的线性规划基本定理的证明。下面要具体说说怎么算。

3.算法步骤

Step1: (当然要先化为标准型)然后找到初始可行基B,解 B x B = b Bx_B=b BxB​=b,求得 x B = B − 1 b = b ~ x_B=B^{-1}b=\tilde{b} xB​=B−1b=b~,令 x N = 0 x_N=0 xN​=0,确定初始可行解。并计算初始目标函数值 f = c B x B , f=c_Bx_B, f=cB​xB​,其中 c B c_B cB​代表 x B x_B xB​对应的目标函数系数向量。

Step2: 求单纯形乘子 w = c B B − 1 w=c_BB^{-1} w=cB​B−1。对于所有非基变量,计算判别(检验)数 ( z j − c j ) = w p j − c j (z_j-c_j)=wp_j-c_j (zj​−cj​)=wpj​−cj​ 。令 z k − c k = m a x j ∈ R N z_k-c_k=max_{j \in R_N} zk​−ck​=maxj∈RN​​{ z j − c j z_j-c_j zj​−cj​}, ( R N (R_N (RN​是非基变量的下标集)若 z k − c k ≤ 0 z_k-c_k \leq 0 zk​−ck​≤0,则对于所有非基变量 z j − c j ≤ 0 z_j-c_j \leq 0 zj​−cj​≤0,而对于基变量的判别数总是零,因此停止计算,现行基本可行解是最优解。否则,进行下一步。

Step3: 解 B k y k = p k B_ky_k=p_k Bk​yk​=pk​,得到 y k = B k − 1 p k y_k=B_k^{-1}p_k yk​=Bk−1​pk​,若 y k ≤ 0 y_k\leq 0 yk​≤0,即 y k y_k yk​的每一个分量均为非正数,则停止计算:问题不存在有限最优解。否则,进行第四步。

Step4: 确定下标 r : b r ~ y r k = m i n r:\frac{\tilde{b_r}}{y_{rk}}=min r:yrk​br​~​​=min{ b i ~ y i k ∣ y i k > 0 \frac{\tilde{b_i}}{y_{ik}}|y_{ik}>0 yik​bi​~​​∣yik​>0}。 x B r x_{Br} xBr​为离基变量, x k x_{k} xk​为进基变量。用 p k p_k pk​替换 p B r p_{Br} pBr​,得到新的矩阵B,返回第一步。

考试计算中,我们经常使用单纯形表完成上述计算。 不懂? 看栗子就完事了! 在这里插入图片描述

4.算例

在这里插入图片描述

ONE: 化为标准型,并找到初始可行基B确定对应的基本可行解 x B x_B xB​ 标准型如下: 在这里插入图片描述

从标准型容易得到一个以 B = I 2 B=I_2 B=I2​为基底的基本可行解: x B = B − 1 b = [ 80 , 90 ] T , x = [ x N , x B ] = [ 0 , 0 , 80 , 90 ] T , c B = [ 0 , 0 ] , x_B=B^{-1}b=[80,90]^T,x=[x_N,x_B]=[0,0,80,90]^T,c_B=[0,0], xB​=B−1b=[80,90]T,x=[xN​,xB​]=[0,0,80,90]T,cB​=[0,0],非基变量的下标集 R N = R_N= RN​={ 1 , 2 1,2 1,2}。

TWO:计算初始检验数 q j q_j qj​及目标函数值 z 0 z_0 z0​,建立初始单纯形表 计算检验数 q j = z j − c j = w p j − c j = c B B − 1 p j − c j q_j=z_j-c_j=wp_j-c_j=c_BB^{-1}p_j-c_j qj​=zj​−cj​=wpj​−cj​=cB​B−1pj​−cj​和目标函数值 z 0 = f = c B T x B = c B B − 1 b z_0=f=c_B^Tx_B=c_BB^{-1}b z0​=f=cBT​xB​=cB​B−1b,建立如下图所示的单纯形表: 在这里插入图片描述 具体对于本题来说这里 w = c B B − 1 = 0 w=c_BB^{-1}=0 w=cB​B−1=0,因此初始检验数为 q j = − c j q_j=-c_j qj​=−cj​,目标函数值为 z 0 = f 0 = c B T x B = 0 , z_0=f_0=c_B^Tx_B=0, z0​=f0​=cBT​xB​=0,因此初始单纯形表如下: 在这里插入图片描述 THREE:决定进基矢量 a k a_k ak​ 取最大检验数 q 2 = 16 > 0 q_2=16>0 q2​=16>0所对应的变量作为进基矢(变)量,即 k = 2 k=2 k=2。

FOUR:决定离基矢量 a r a_r ar​和主元素 y r k y_{rk} yrk​ y i j , j ≠ 0 y_{ij},j\not=0 yij​,j​=0即单纯形表第i行的第j个变量对应列的对应元素, y i 0 y_{i0} yi0​即单纯形表最后一列( B − 1 b B^{-1}b B−1b)的第i行元素。 首先计算比值 y i 0 / y i k y_{i0}/y_ik yi0​/yi​k,这里 k = 2 , i = 1 , 2 k=2,i=1,2 k=2,i=1,2。结果为: y 10 y 12 = 80 / 4 = 20 ; \frac{y_{10}}{y_{12}}=80/4=20; y12​y10​​=80/4=20; y 20 y 22 = 90 / 3 = 30 \frac{y_{20}}{y_{22}}=90/3=30 y22​y20​​=90/3=30 故由 r : b r ~ y r k = m i n r:\frac{\tilde{b_r}}{y_{rk}}=min r:yrk​br​~​​=min{ b i ~ y i k ∣ y i k > 0 \frac{\tilde{b_i}}{y_{ik}}|y_{ik}>0 yik​bi​~​​∣yik​>0}。得 r = 1 r=1 r=1。因此主元素为 y 12 y_{12} y12​,离基变量为 x 3 x_3 x3​。

FIVE:以 y r k y_{rk} yrk​为主元素更新(进行Gauss消元)单纯形表 对单纯形表进行初等行变换使得主元素 y r k y_{rk} yrk​所在位置变为1,且对应列的其他行都变为0。 对于本题将第一行除以4,然后分别乘以(-3)和(-16)加到第2、3行,完成第一次迭代,函数值减小为为-320。 在这里插入图片描述 然后回到THREE进行下次迭代。 那么什么时候停止呢?下面来说 在这里插入图片描述

5.算法收敛性

以极小化问题为例,最大检验数 q r = m a x j ∈ R N ( z j − c j ) = m a x j ∈ R N ( w p j − c j ) = m a x j ∈ R N ( c B B − 1 p j − c j ) q_r=max_{j\in R^N}(z_j-c_j)=max_{j\in R^N}(wp_j-c_j)=max_{j\in R^N}(c_BB^{-1}p_j-c_j) qr​=maxj∈RN​(zj​−cj​)=maxj∈RN​(wpj​−cj​)=maxj∈RN​(cB​B−1pj​−cj​)每次迭代必然出现下面三种情形: (1) q r ≤ 0 q_r\leq 0 qr​≤0:这时现行基本可行解就是最优解。 (2) q r > 0 , y r k ≤ 0 q_r> 0,y_{rk}\leq0 qr​>0,yrk​≤0:当 x k x_k xk​无线增大时,目标函数趋于负无穷,因此解无解。 (3) q r > 0 , y r k > 0 q_r> 0,y_{rk}>0 qr​>0,yrk​>0:这时求出的新的基本可行解,经迭代使目标函数下降。 收敛定理: 对于非退化问题(即存在非退化基本可行解),单纯形方法经有限次迭代后达到最优基本可行解,或得出无界的结论。 所以上面的那个例子在迭代了2次后最大检验数 q r ≤ 0 q_r\leq 0 qr​≤0,迭代停止,此时的解即为最优解。 在这里插入图片描述

6.Matlab实现

这里实现已经化为标准型且约束方程是规范形式。其实后面学习我们知道其他线性优化问题均可以通过M法或两阶段法转化至此。 通过上述算例的步骤一步一步实现即可,不算困难。大致可分四步:1.输入问题;2.建立初始单纯形表;3.迭代寻找最优解;4.输出结果。 因为这次期末时间过于紧张,这里先把代码总表给出,之后再按下面的步骤走一遍。 在这里插入图片描述

1.输入问题 2.建立初始单纯形表 3.迭代寻找最优解 4.输出结果 5.附录(代码总表) % 输入一标准型,实现程序求求解 % 若有m个方程,则初始基本可行解为E,且位于最后mxm矩阵 % matlab代码 % Made by Fei % 返回最优解(变量名和对应取值)和最优函数值 %c,A,b分别为目标函数系数向量,约束函数系数矩阵,约束函数右端常数 function [Y,T] = Mysolve(c,A,b) % 预备工作 % 求出变量个数以及方程个数 n_numx = size(A,2); n_eq = rank(A); % 记录初始基本可行解及相应的变量所在位置 XB = b; XB_index=[]; for i=1:n_eq XB_index(i)=n_numx-n_eq+i; end % 求出单纯形乘子 cB = c(:,(n_numx-n_eq+1):n_numx); w = cB; % 求出初始函数值 f0 = cB*XB; % 求出n_numx个检验数 z=[]; for i=1:n_numx if(i>n_numx-n_eq) z(i)=0; else temp = w*A(:,(n_numx-n_eq-i+1))-c(:,i); z(i)=temp; end end % 开始迭代求解 flag=0; while(flag==0) % 判断检验数是否全部全部为非正数,如是则停止迭代平输出结果 flag=1; for i=1:n_numx temp=z(i); if(temp>0) flag=0; end end % 如果全负则输出结果 if(flag==1) Y=[XB_index;XB']; % 变量名及所对应的取值 T=f0; % 最优目标函数值 break end %求解过程 % 找出最大检验数及并确定进基变量 Max=z(1); index_in=1; for i=2:n_numx if(z(i)>Max) Max=z(i); index_in=i; end end % 找出yr以及出基变量 y=b./A(:,index_in); index_out_H=0; for i=1:n_eq if(y(i)


【本文地址】

公司简介

联系我们

今日新闻

    推荐新闻

    专题文章
      CopyRight 2018-2019 实验室设备网 版权所有