单纯形法求解线性规划

Galerkin码农选手

已于 2022-08-03 08:52:16 修改

阅读量1k

点赞数

分类专栏：优化方法文章标签：动态规划

于 2021-06-26 09:13:10 首次发布

本文链接：https://blog.csdn.net/forrestguang/article/details/118113065

版权

优化方法专栏收录该内容

12 篇文章 24 订阅

订阅专栏

说明：以下内容为本人上机作业，仅供读者学习参考使用，禁止抄袭

问题描述

考虑从源节点到端节点传输数据的网络，我们首先引入一些定义：
$\bullet$ K\quad flows：传输K中不同的数据。
$\bullet$ L\quad link：连接网络节点的L条无向边。
$\bullet$ $c_l$ ：在第 $l$ 条边上传输数据的最大容量。
$\bullet$ $P_k$ ：对于不同的k flow，运输成功可选择的路path的数目。
$\bullet$ $\mathbf{x_k}=\left(x_{k,1}, \ldots, x_{k,P_{K}}\right)^T$ ： $x_{k,i}$ 表示对于k\quad flow，在第 $i$ 条路某条边分配的数据量。
给定向量 $C=\left(c_{1}, c_{2}, \ldots, c_{L}\right)^T$ ，其中 $c_l \in R^{+},\forall{l}$ 。定义：
$R_{k}=\left(\begin{array}{cccc} R_{1,1}^{k} & R_{1,2}^{k} & \cdots & R_{1, P_{k}}^{k} \\ R_{2,1}^{k} & R_{2,2}^{k} & \cdots & R_{2, P_{k}}^{k} \\ \vdots & \vdots & \ddots & \vdots \\ R_{L, 1}^{k} & R_{L, 2}^{k} & \cdots & R_{L, P_{k}}^{k} \end{array}\right)$
其中 $R_{k}$ 中的任意一个元素 $\in \{0,1\}$ ， $P=\sum_{i=1}^{K}P_i$ ，定义 $\in R^{L\times P}$ ，其中 $R_{l,i}^k = 1$ 表示运输k\quad flow时，第 $i$ 条路会经过第 $l$ 条边：
$R=\left(R_{1}, \ldots, R_{K}\right)$
定义向量 $x_k=\left(x_{k,1}, x_{k,2}, \ldots, x_{k,P_{k}}\right)^T$ ，其中 $x_{k,i} \geq 0,\forall{k,i}$ ,且 $x_{k,P_{k}} > 0,\forall{k}$ 。定义 $(x_1,\ldots,x_k)\in R^P$ 。有了这些定义以后，给出Problem 1的描述：
$\min _{X} \max _{l} \frac{R[l] x}{c_{l}}$
$\text { s.t. } Rx \leq C \\{ } \left\|x_{k}\right\|_{1}=d_{k} \\ x \geq 0$
转化为标准线性规划问题
引入中间变量 $\max_{l}\frac{y_l}{c_l}$ ，则可以写成下面形式:
$\min _{t,x,\alpha,\beta} t$
$\text { s.t. } Rx + \alpha = C,\\{ } Rx + \beta = tC, \\{ } \boldsymbol{1}^T x_k = d_k,\\ x \geq 0,\alpha \geq 0,\beta \geq 0,t \geq 0.$
进一步修改得到 $\longrightarrow$ ：
$\min _{t,x,\alpha,\beta} t,$
$\text { s.t. } Rx + \alpha = C,\\{ } tC + \alpha - \beta = C, \\{ } \boldsymbol{1}^T x_k = d_k,\\ x \geq 0,\alpha \geq 0,\beta \geq 0,t \geq 0.$
这样我们引入定义:\
$\bullet$ $\hat{c} = (0,\ldots,0,1)$ \
$\bullet$ $\hat{x} = (x,\alpha,\beta,t)^T$ \
$\bullet$ $(c_1,\ldots,c_l,c_L,\ldots,c_L,d_1,\ldots,d_K)^T$ \
得到最终的标准线性规划问题：
$\min _{\hat{x}} \hat{c}^T\hat{x},$
$\text { s.t. } A \hat{x} = b,\\{ } \hat{x} \geq 0.$
其中定义 $I_L$ 为 $L$ 阶单位矩阵， $(c_1,\ldots,c_L)^T$ 以及 $S\in R^{K\times P}$ ，矩阵 $S$ 的第一行只有前 $P_1$ 个元素为1，其余为0，第二行元素只有 $P_1 + 1,\ldots,P_1 + P_2$ 为1，其余为0，类似地，最后一行只有后面 $P_K$ 个元素为1，其余为0.由此得到下面这个系数矩阵 $A$ 的表达式：
$A=\left[\begin{array}{cccc} R & I_L & 0 & 0 \\ 0 & I_L & -I_L & C \\ S & 0 & 0 & 0 \end{array}\right]$

标准线性规划

考虑标准的线性规划问题，其中 $\in R^{m \times n},C \in R^{n},b \in R^{m}$ ：
$min_{x} C^{T}x$
$\text { s.t. } Ax = b,$
$\geq 0.$
众所周知，在利用单纯形法求解线性规划的时候，从当前极点 $x$ 到下一个极点 $x^{+}$ 的迭代中，假设 $x$ 的基指标为 $B$ ，非基指标为 $N$ ，我们有
$c^Tx^{+}=c^Tx + (C_N - N^{T}B^{-T}C_B)^Tx_{N}^{+}$
由于在指标 $N$ 中，根据单纯形法， $x^{+}$ 只有一个元素大于0，记为 $x_{q}^{+}$ ，如果
$C_N - N^{T}B^{-T}C_B)^Tx_{N}^{+} < 0$
则 $x$ 非最优点，再根据
$x_{B}^{+}=B^{-1}b-B^{-1}A_qx_{q}^{+}$
我们需要选取适当的 $x_{q}^{+}$ ，由此引出了以下几种准则。
Dantzig准则
$min_{j \in N}(C_N - N^{T}B^{-T}C_B)_j,(C_N - N^{T}B^{-T}C_B)_j < 0.$
关于这个Dantzig准则，我们只选取了具体的需要修改的索引位置 $q$ ，但是具体的 $x_{q}^{+}$ 取值为多少，还需要根据来确定，具体的做法是选尽量大的 $x_{q}^{+}$ ，使得 $x_{B}^{+}$ 只有一个元素为0，其他元素均大于或等于0。
Steepest Edge 准则
$min_{j \in N}\frac{(C_N - N^{T}B^{-T}C_B)_j}{B^{-1}N_j},(C_N - N^{T}B^{-T}C_B)_j < 0.$
我们可以这么理解：上面的Dantzig准则类似于最速下降法，每次迭代选取当前迭代的最优方向，但是整体来看却不见得最优，这个Steepest Edge 准则有点类似于做了一个正规化。

大M法求解初始基本可行解

在求解线性规划问题时，首先我们需要一个初始点满足约束，然后从初始点出发进行单纯形法迭代，经过调研，大M法和两阶段法比较普遍，下面我们重点介绍大M法求初始基本可行解的过程.
我们考虑标准的线性规划问题，引入一个充分大的常数M，定义 $\hat{C} = (c_1,\ldots,c_n,M,\ldots,M),\hat{A} = (A,I_{m}),\hat{x} = (x,x_{n + 1},\ldots,x_{n + m})$ ，现在考虑问题：
$\min _{X} \hat{C}^{T}\hat{x}$
$\text { s.t. } \hat{A}\hat{x} = b,$
$\hat{x} \geq 0.$
则对于线性规划问题，我们就有第一个初始可行点 $\hat{x}=(0,\ldots,0,b1,\ldots,b_m)$ ，有了这个初始可行点以后，我们就可以根据单纯形法做迭代，得到一个形式为 $(x_1,\ldots,x_n,0,\ldots,0)$ 的解，最终得到一个原始问题的初始可行点 $(x_1,\ldots,x_n)$ 。
说明：这种方法仅限于约束项的 $\geq 0$ 的情况(这样才保证大M法在第一步是对的)，这里我们需要把M选取为充分大的常数，这样就可以保证，如果原始问题有可行点，那么新问题的最优解一定是类似于 $(x_1,\ldots,x_n,0,\ldots,0)$ 的形式，在我们具体操作过程中，这个M的选取并不固定，本人做了几个实验，发现如果采取高斯分布生成 $A, b, c$ ，由于随机种子的不一样，有的时候M的选取会超过 $10^{20}$ ，不过在我们这次项目中， $A, b, c$ 的元素相对较小，所以不会受到影响。
一般线性规划大M法的代码，求解初始可行解

import numpy as np
import time
import matplotlib.pyplot as plt
def Bland(s,indn):
    ind = np.where(s[:,0]<0)[0]
    q = ind[0]
    return indn[q], q

def DanzigSelect(s,indn):#indn非基索引
    q = s.argmin()
    return indn[q], q

def SteepSelect(A_n,s,indn):
    ratio = s/A_n.transpose()
    q = ratio.argmin()
    return indn[q], q
def bigM(c,A,b,piv,M):#使用大M方法求解初始基本可行解
    eps = 1e-6
    [m,n] = A.shape
    x = np.zeros([n + m,1]);x[n:,:] = b
    
    Ah = np.hstack([A,np.identity(m)])
    ch =  np.vstack([c,M*np.ones([m,1])])
    indb = np.arange(m + n)[n:]#初始基指标索引
    index = np.arange(0,m + n)
    indn = np.delete(index,indb)
    B = Ah[:,indb]#初始B是单位矩阵
    B_inv = B
    #print(B)
    N = Ah[:,indn]
    iteration = 0        # determine the max iteration
    
    while iteration < 500:
        y = np.dot(B_inv.T, ch[indb,:])
        s = ch[indn,:] - N.T@y#C_N - N^TB^{-T}C_B
        if all(s >= 0):#说明这个就是最优解
            print('unbelievabel,the resid:%.2e, done,the epoch:%d'%(np.linalg.norm(Ah@x - b),iteration))
            break
        if piv == 'Danzig':
            q1, q2 = DanzigSelect(s, indn)
        if piv == 'steepest-edge':
            A_re = B_inv@N
            A_n = np.linalg.norm(A_re,axis = 0,keepdims = True)
            q1, q2 = SteepSelect(A_n, s, indn)
        if piv == 'Bland':
            q1, q2 = Bland(s, indn)
        #q1是非基指标中需要被替换的索引
        u = B_inv@Ah[:, q1:(q1+1)]   
        if all(u <= 0):#说明最优值无界
            print('Linear program is unbounded.')
            break
            
        inds1 = np.where(u[:,0] > 0)[0]#找出所有大于0的索引
        inds2 = indb[inds1]#在基指标中把对应大于0的索引全找出来
        ratio = x[inds2,:]/u[inds1,:]
        x_q = min(ratio)        # a number
        p1 = inds2[ratio.argmin()]
        p2 = inds1[ratio.argmin()]
        
        x[indb,:] = x[indb,:]- x_q*u
        vec = np.zeros([indn.size,1])
        vec[q2,0] = x_q
        x[indn,:] = x[indn,:] + vec
        rho = (x[n:,:]**2).mean()
        if iteration%50 == 0:
            
            print('epoch:%d,enter basis:%d,quit basis:%d,resid:%.2e,rest err:%.2e'
                  %(iteration,q1,p1,np.linalg.norm(Ah@x - b),rho))
        
        if rho < eps:
            print('the end epoch:%d,rho:%.2e,the real resid:%.2e'
                  %(iteration,rho,np.linalg.norm(A@x[:n,:] - b)))
            break
        else:
            indb = np.delete(indb,p2)
            indb = np.append(indb,q1)
            indn = np.delete(indn,q2)
            indn = np.append(indn,p1)
            
            B = Ah[:, indb]
            B_inv = np.linalg.inv(B)
            N = Ah[:, indn]
            iteration = iteration + 1
    return x[:n,:]

def simplex(c,A,b,piv,M):
    eps = 1e-6
    [m,n] = A.shape
    x = bigM(c,A,b,piv,M)
    #print(x.shape,A.shape,b.shape,indb.size,m)
    if np.linalg.norm(A@x - b) > eps or (any(x < 0) and min(x) < -eps):
        print('the init solution error')
        return 0
    else:
        print('the init solution come on')
    indb = []
    geq = np.where(x > 0)[0]
    eq = np.where(x <= 0)[0]
    indb = np.append(indb,geq)
    if indb.size == m:
        pass
    elif indb.size < m:
        er = m - indb.size
        print('basic need:%d'%(er))
        ra = np.random.choice(eq,er,replace = False)
        indb = np.append(indb,ra)
    else:
        print('wrong')
        return 0
    index = np.arange(0,n)
    indb = indb.astype('int32')
    indn = np.delete(index,indb)
    B = A[:,indb]
    
    B_inv = np.linalg.inv(B)
    
    N = A[:,indn]
    iteration = 0        # determine the max iteration
    
    while iteration < 2000:
        y = np.dot(B_inv.T, c[indb,:])
        s = c[indn,:] - N.T@y#C_N - N^TB^{-T}C_B
        value = (c*x).sum()
        if all(s >= 0) or (any(s < 0) and min(s) > -eps):#说明这个就是最优解
            print('the end epoch:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,np.linalg.norm(A@x - b),value))
            break
        if piv == 'Danzig':
            q1, q2 = DanzigSelect(s, indn)
        if piv == 'steepest-edge':
            A_re = B_inv@N
            A_n = np.linalg.norm(A_re,axis = 0,keepdims = True)
            q1, q2 = SteepSelect(A_n, s, indn)
        if piv == 'Bland':
            q1, q2 = Bland(s, indn)
        #q1是非基指标中需要被替换的索引
        u = B_inv@A[:, q1:(q1+1)]   
        if all(u <= 0):#说明最优值无界
            print('Linear program is unbounded.')
            break
            
        inds1 = np.where(u[:,0] > 0)[0]#找出所有大于0的索引
        inds2 = indb[inds1]#在基指标中把对应大于0的索引全找出来
        ratio = x[inds2,:]/u[inds1,:]
        x_q = min(ratio)        # a number
        p1 = inds2[ratio.argmin()]
        p2 = inds1[ratio.argmin()]
        
        x[indb,:] = x[indb,:]- x_q*u
        vec = np.zeros([indn.size,1])
        vec[q2,0] = x_q
        x[indn,:] = x[indn,:] + vec
        res = np.linalg.norm(A@x - b)
        if res > eps:
            print('blow up,epoch:%d,resid:%.2e'%(iteration,res))
            break
        if iteration%100 == 0:
            #print(np.linalg.det(B))
            print('epoch:%d,enter basis:%d,quit basis:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,q1,p1,res,value))
        indb = np.delete(indb,p2)
        indb = np.append(indb,q1)
            
        indn = np.delete(indn,q2)
        indn = np.append(indn,p1)
            
        B = A[:, indb]
        B_inv = invB(B_inv,p2,u)
           
        N = A[:, indn]
        iteration = iteration + 1
    return x

L1 = 100
P1 = 200
k = 5
p = 0.02

np.random.seed(1234)

R = np.random.binomial(1,p,[L1,P1])#二项分布，生成[L1,P1]行数组，元素全是0或1
d = np.random.gamma(2,2,[k,1])*1#gamma分布，
u = np.random.gamma(2,2,[L1,1])*3+20
s = P1 - k
ind = np.random.randint(0,k,s)#随机抽取0到K之间的整数
S = np.zeros((k,s))
for i in range(s):
    j = ind[i]
    S[j,i] = 1
S = np.hstack((np.identity(k),S))

n = P1+2*L1+1       #364
m = 2*L1+k          #196

# construct the vector c, only the last element is 1
c = np.zeros((n,1))
c[-1,0] = 1

print(c.shape, c[-1,0])

# construct the matrix A
S1 = np.hstack((R,np.zeros((L1,L1)),np.identity(L1),np.zeros((L1,1))))
S2 = np.hstack((np.zeros((L1,P1)),-np.identity(L1),np.identity(L1),u))
S3 = np.hstack((S,np.zeros((k,(2*L1+1)))))
A = np.vstack((S1,S2,S3))

b = np.vstack((u,u,d))
piv = 'Bland'

M = 1e6
x = bigM(c,A,b,piv,M)
#x = simplex(c,A,b,piv,M)

在这里插入图片描述

两阶段法求解初始可行解

简单介绍一下两阶段法的做法，和上面提到的大M法相比，两阶段法最大的优势就是不需要选择超参数M。仍然以标准的线性规划问题举例，我们的求解过程分成两个阶段进行。下面考虑第一个阶段，其中 $(y_1,\ldots,y_m)$ ：
$\min_{x,y} \sum_{i = 1}^m y_i$
$\text { s.t. } Ax + I_m y = b,$
$\geq 0.$
求解线性规划，如果求出最优值不为0，说明原始线性规划没有可行解。如果最优值求出来为0，那么我们需要分情况考虑：
第一种情况：问题的最优解 $(x, y)$ 的基指标包含人工变量 $y$ 的指标，此时说明原始线性规划退化，这个时候需要做一些替换，不断把包含 $y$ 的基指标替换出来，然后得到初始可行解。
第二种情况：问题的最优解 $(x, y)$ 的基指标不包含人工变量 $y$ 的指标，此时说明我们可以直接把 $x$ 作为一个初始可行解来使用。
在第一个阶段的基础上进行第二个阶段就很简单，直接根据初始可行解做单纯形法迭代即可。
一般线性规划两阶段法代码

import numpy as np

def Bland(s,indn):
    ind = np.where(s[:,0]<0)[0]
    q = ind[0]
    return indn[q], q

def DanzigSelect(s,indn):#indn非基索引
    q = s.argmin()
    return indn[q], q

def SteepSelect(A_n,s,indn):
    ratio = s/A_n.transpose()
    q = ratio.argmin()
    return indn[q], q
def invB(B_inv,p,u):
    Bv = B_inv.copy()
    Bv[p,:] = Bv[p,:]/u[p,0]
    m = u.size
    ind = []
    for i in range(m):
        if i == p:
            pass
        else:
            ind.append(i)
            Bv[i,:] = Bv[i,:] - u[i,0]*Bv[p,:]
    ind.append(p)
    Bv = Bv[ind,:]
    return Bv
def stage_one(A,b,piv):#使用大M方法求解初始基本可行解
    eps = 1e-6
    [m,n] = A.shape
    x = np.zeros([n + m,1]);x[n:,:] = b
    
    Ah = np.hstack([A,np.identity(m)])
    ch =  np.zeros([n + m,1]);ch[n:,:] = np.ones([m,1])
    indb = np.arange(m + n)[n:]#初始基指标索引
    index = np.arange(0,m + n)
    indn = np.delete(index,indb)
    B = Ah[:,indb]#初始B是单位矩阵
    B_inv = B
    
    N = Ah[:,indn]
    iteration = 0        # determine the max iteration
    
    while iteration < 1000:
        
        y = np.dot(B_inv.T, ch[indb,:])
        s = ch[indn,:] - N.T@y#C_N - N^TB^{-T}C_B
        if all(s >= 0):#说明这个就是最优解
            print('unbelievabel,the resid:%.2e, done,the epoch:%d'%(np.linalg.norm(Ah@x - b),iteration))
            break
        if piv == 'Danzig':
            q1, q2 = DanzigSelect(s, indn)
        if piv == 'steepest-edge':
            A_re = B_inv@N
            A_n = np.linalg.norm(A_re,axis = 0,keepdims = True)
            q1, q2 = SteepSelect(A_n, s, indn)
        if piv == 'Bland':
            q1, q2 = Bland(s, indn)
        #q1是非基指标中需要被替换的索引
        u = B_inv@Ah[:, q1:(q1+1)]   
        if all(u <= 0):#说明最优值无界
            print('Linear program is unbounded.')
            break
            
        inds1 = np.where(u[:,0] > 0)[0]#找出所有大于0的索引
        inds2 = indb[inds1]#在基指标中把对应大于0的索引全找出来
        ratio = x[inds2,:]/u[inds1,:]
        x_q = min(ratio)        # a number
        p1 = inds2[ratio.argmin()]
        p2 = inds1[ratio.argmin()]
        
        x[indb,:] = x[indb,:]- x_q*u
        vec = np.zeros([indn.size,1])
        vec[q2,0] = x_q
        x[indn,:] = x[indn,:] + vec
        aug = np.linalg.norm(Ah@x - b)
        if aug > eps:
            print('blow up,epoch:%d,aug resid:%.2e'%(iteration,aug))
            break
        value = x[n:,:].sum()
        if iteration%100 == 0:
            
            print('epoch:%d,enter basis:%d,quit basis:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,q1,p1,aug,value))
        
        if value < eps:
            print('the end epoch:%d,opt val:%.2e,the real resid:%.2e'
                  %(iteration,value,np.linalg.norm(A@x[:n,:] - b)))
            break
        else:
            
            indb = np.delete(indb,p2)
            indb = np.append(indb,q1)
            
            indn = np.delete(indn,q2)
            indn = np.append(indn,p1)
            
            B = Ah[:, indb]
            B_inv = invB(B_inv,p2,u)
           
            N = Ah[:, indn]
            iteration = iteration + 1
    return x[:n,:]

def stage_two(c,A,b,piv,L,P,K):
    eps = 1e-6
    [m,n] = A.shape
    x = stage_one(A,b,piv)
    
    if np.linalg.norm(A@x - b) > eps or (any(x < 0) and min(x) < -eps):
        print('the init solution error')
        print('---------------------------------------')
        return 0
    else:
        print('the init solution come on')
        print('---------------------------------------')
    geq = np.where(x > eps)[0]
    eq = np.where(x == 0)[0]
    #eeq = np.where(any(x > 0 and x < eps))[0]
    if geq.size < m:
        er = m - geq.size
        ra = np.random.choice(eq.size,er)
        indb = np.append(geq,eq[ra])
    elif geq.size > m:
        print('error')
        return 0
    else:
        indb = geq
    index = np.arange(0,n)
    indn = np.delete(index,indb)
    B = A[:,indb]
    
    B_inv = np.linalg.inv(B)
    
    N = A[:,indn]
    iteration = 0        # determine the max iteration
    
    while iteration < 2000:
        y = np.dot(B_inv.T, c[indb,:])
        s = c[indn,:] - N.T@y#C_N - N^TB^{-T}C_B
        value = (c*x).sum()
        if all(s >= 0) or (any(s < 0) and min(s) > -eps):#说明这个就是最优解
            print('the end epoch:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,np.linalg.norm(A@x - b),value))
            break
        if piv == 'Danzig':
            q1, q2 = DanzigSelect(s, indn)
        if piv == 'steepest-edge':
            A_re = B_inv@N
            A_n = np.linalg.norm(A_re,axis = 0,keepdims = True)
            q1, q2 = SteepSelect(A_n, s, indn)
        if piv == 'Bland':
            q1, q2 = Bland(s, indn)
        #q1是非基指标中需要被替换的索引
        u = B_inv@A[:, q1:(q1+1)]   
        if all(u <= 0):#说明最优值无界
            print('Linear program is unbounded.')
            break
            
        inds1 = np.where(u[:,0] > 0)[0]#找出所有大于0的索引
        inds2 = indb[inds1]#在基指标中把对应大于0的索引全找出来
        ratio = x[inds2,:]/u[inds1,:]
        x_q = min(ratio)        # a number
        p1 = inds2[ratio.argmin()]
        p2 = inds1[ratio.argmin()]
        
        x[indb,:] = x[indb,:]- x_q*u
        vec = np.zeros([indn.size,1])
        vec[q2,0] = x_q
        x[indn,:] = x[indn,:] + vec
        res = np.linalg.norm(A@x - b)
        if res > eps:
            print('blow up,epoch:%d,resid:%.2e'%(iteration,res))
            break
        if iteration%100 == 0:
            #print(np.linalg.det(B))
            print('epoch:%d,enter basis:%d,quit basis:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,q1,p1,res,value))
        indb = np.delete(indb,p2)
        indb = np.append(indb,q1)
            
        indn = np.delete(indn,q2)
        indn = np.append(indn,p1)
            
        B = A[:, indb]
        B_inv = invB(B_inv,p2,u)
           
        N = A[:, indn]
        iteration = iteration + 1
    return x


L = 64
P = 128
K = 30
p = 0.02

np.random.seed(1)

# seed 1 with 512,1024, Danzig beats steepest
# seed 1 with 256, 512 steepest beats Danzig

R = np.random.binomial(1,p,[L,P])#二项分布，生成[L1,P1]行数组，元素全是0或1
d = np.random.gamma(2,2,[K,1])*1#gamma分布，
u = np.random.gamma(2,2,[L,1])*3+20
s = P - K
ind = np.random.randint(0,K,s)#随机抽取0到K之间的整数
S = np.zeros((K,s))
for i in range(s):
    j = ind[i]
    S[j,i] = 1
S = np.hstack((np.identity(K),S))

#---------------------------------
L = R.shape[0];P = R.shape[1];K = S.shape[0]
m = 2*L + K;n = P + 2*L + 1
print(type(R),type(S),type(u),type(d))
print(R.shape,S.shape,u.shape,d.shape)
print(max(R.reshape(-1,1)),min(R.reshape(-1,1)),max(S.reshape(-1,1)),min(S.reshape(-1,1)))
row1 = np.hstack([R,np.identity(L),np.zeros([L,L]),np.zeros([L,1])])
row2 = np.hstack([np.zeros([L,P]),np.identity(L),-np.identity(L),u])
row3 = np.hstack([S,np.zeros([K,2*L + 1])])
A = np.vstack([row1,row2,row3])
print(np.linalg.matrix_rank(A),m)
b = np.vstack([u,u,d])
c = np.zeros([n,1]);c[-1,0] = 1.0
#piv = 'Bland'
piv = 'Danzig'
x = stage_two(c,A,b,piv,L,P,K)


print(all(x >= 0))
geq = np.where(x > 0)[0]
eq = np.where(x == 0)[0]
leq = np.where(x < 0)[0]
print(geq.size,eq.size,leq.size,m)
print('the nagetive value:')
print(x[leq])
#----------------------

在这里插入图片描述

单纯形法求解上述线性规划

初始可行解的选取
在我们接连介绍了单纯形法的一些准则，概念，以及给出了我们项目需要求解的线性规划具体的形式以后，我们将在这一小节给出使用单纯形法求解我们的项目的具体做法。
首先，我们需要获得第一个初始可行解来启动单纯形法，在实际操作种，我们发现如果简单利用大M法或者是两阶段法来求解第一个初始可行解，当矩阵规模过大的时候，算法变得很不稳定，由此我们想出了一种专门针对这样形式的矩阵的单纯形构造法来构造第一个可行解。具体做法是，首先给出第一个子问题如下：
$\begin{array}{cl} \min _{x,y} & \sum_{i = 1}^{K + L} z_i， \\ \text { s.t. } & Sx = d,\\{ } & Rx + y = C,\\{ }& x \geq 0. \end{array}$
利用两阶段法的第一个阶段来求解第一个子问题\ref{sub1}，不妨记为 $x,y)=(x_{1:P},x_{P:P + L})$ 。\
然后我们定义一个变量 $t$ ，定义后续给出。现在我们定义第二个子问题： $x_{P:P+L} - x + tC=C$
求解第二个子问题我们得到的解记为 $x_{P+L:P+2L}$ ，很容易看出 $x_{P+L:P+2L}=x_{P:P+L} + (t-1)C$ 。最终我们得到第一个初始可行解为：
$\bullet$ $x = (x_{1:P},x_{P:P+L},x_{P+L:P+2L},t)$
根据上面求解过程，我们知道 $x_{1:P} \geq 0，x_{P:P+L} \geq 0$ ，并且元素为0 的个数至少是 $P - K$ 。如果元素为0的个数是 $P - K + 1$ ,我们选取合适的 $\geq 0$ ，则会有 $\geq 0$ ，此时得到的解非0元素的个数至少是 $P + 1 - K = m - n$ ，那么这个就可以作为我们的初始可行解。\par
如果元素为0的个数恰好是 $P - K$ ，我们可以选取合适的 $\geq 0$ ，这里选择 $\min_{0 \leq j \leq L-1}\frac{x_{P:P+j}}{C_j}$ ，这个其实还原为原始问题就是 $\min_{0 \leq j \leq L-1}\frac{\alpha}{C_j}$ ，由于 $Rx+\alpha=C$ ，很明显有 $\geq 0$ ，使得 $x_{P+L:P+2L}$ 恰好只有一个元素是0，其余都大于0，经过这个过程，得到的解非0元素的个数至少是 $P + 1 - K = n - m$ ，那么也得到初始可行解
迭代法的选择
在得到初始可行解以后，我们需要做一个判断，如果初始可行解 $x_0$ 存在 $m$ 个元素大于0，其余元素为0，那么我们就可以得到相应的基指标。如果初始可行解 $x_0$ 大于0的元素个数小于 $m$ ，那么我们需要把那些大于0元素的索引取出来，然后再根据适当的原则选取差额索引数目放入得到基指标。
在具体的代码实现中，我们可以有上述提到的两种准则可以使用，这里需要求解基矩阵 $B$ 的逆矩阵，这个我们也做了相应的处理。在给定上一次迭代使用的基矩阵我们记为 $B_1$ ，假设 $B_1 = (A_{1},\ldots,A_{p-1},A_{p},A_{p+1},\ldots,A_{m})$ 和对应的逆矩阵 $B_1^{-1}$ 前提下，新的基矩阵 $B_2 = (A_{1},\ldots,A_{p-1},A_{p+1},\ldots,A_{m},A_q)$ ，则我们有 $B_1^{-1}B_2 = (e_1,\ldots,e_{p-1},e_{p+1},\ldots,e_m,B_1^{-1}A_q)$ ，此时我们可以在 $B_1^{-1}$ 的基础上做初等变换得到新矩阵 $B_2$ 的逆矩阵。

数据生成

为了方便展示结果，这里我们展示数据的生成模式：
在这里插入图片描述
上述特定线性规划问题求解的代码
randsimplex.py
随机选择准则更新基指标

import numpy as np
import time
def Bland(s,indn):
    ind = np.where(s[:,0]<0)[0]
    q = ind[0]
    return indn[q], q

def DanzigSelect(s,indn):#indn非基索引
    q = s.argmin()
    return indn[q], q

def SteepSelect(A_n,s,indn):
    ratio = s/A_n.transpose()
    q = ratio.argmin()
    return indn[q], q
def invB(B_inv,p,u):
    Bv = B_inv.copy()
    Bv[p,:] = Bv[p,:]/u[p,0]
    m = u.size
    ind = []
    for i in range(m):
        if i == p:
            pass
        else:
            ind.append(i)
            Bv[i,:] = Bv[i,:] - u[i,0]*Bv[p,:]
    ind.append(p)
    Bv = Bv[ind,:]
    return Bv
def stage_one(A,b):#使用大M方法求解初始基本可行解
    eps = 1e-6
    [m,n] = A.shape
    x = np.zeros([n + m,1]);x[n:,:] = b
    
    Ah = np.hstack([A,np.identity(m)])
    ch =  np.zeros([n + m,1]);ch[n:,:] = np.ones([m,1])
    indb = np.arange(m + n)[n:]#初始基指标索引
    index = np.arange(0,m + n)
    indn = np.delete(index,indb)
    B = Ah[:,indb]#初始B是单位矩阵
    B_inv = B
    
    N = Ah[:,indn]
    iteration = 0        # determine the max iteration
    
    while iteration < 1000:
        
        y = np.dot(B_inv.T, ch[indb,:])
        s = ch[indn,:] - N.T@y#C_N - N^TB^{-T}C_B
        if all(s >= 0):#说明这个就是最优解
            print('unbelievabel,the resid:%.2e, done,the epoch:%d'%(np.linalg.norm(Ah@x - b),iteration))
            break
        z = np.random.rand(1)
        if (z >= 0 and z <= 1/3):
            q1, q2 = DanzigSelect(s, indn)
        elif (z >= 1/3 and z <= 2/3):
            A_re = B_inv@N
            A_n = np.linalg.norm(A_re,axis = 0,keepdims = True)
            q1, q2 = SteepSelect(A_n, s, indn)
        else:
            q1, q2 = Bland(s, indn)
        #q1是非基指标中需要被替换的索引
        u = B_inv@Ah[:, q1:(q1+1)]   
        if all(u <= 0):#说明最优值无界
            print('Linear program is unbounded.')
            break
            
        inds1 = np.where(u[:,0] > 0)[0]#找出所有大于0的索引
        inds2 = indb[inds1]#在基指标中把对应大于0的索引全找出来
        ratio = x[inds2,:]/u[inds1,:]
        x_q = min(ratio)        # a number
        p1 = inds2[ratio.argmin()]
        p2 = inds1[ratio.argmin()]
        
        x[indb,:] = x[indb,:]- x_q*u
        vec = np.zeros([indn.size,1])
        vec[q2,0] = x_q
        x[indn,:] = x[indn,:] + vec
        aug = np.linalg.norm(Ah@x - b)
        if aug > eps:
            print('blow up,epoch:%d,aug resid:%.2e'%(iteration,aug))
            break
        value = x[n:,:].sum()
        if iteration%100 == 0:
            
            print('epoch:%d,enter basis:%d,quit basis:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,q1,p1,aug,value))
        
        if value < eps and all(s) >= 0:
            print('the end epoch:%d,opt val:%.2e,the real resid:%.2e'
                  %(iteration,value,np.linalg.norm(A@x[:n,:] - b)))
            break
        else:
            
            indb = np.delete(indb,p2)
            indb = np.append(indb,q1)
            
            indn = np.delete(indn,q2)
            indn = np.append(indn,p1)
            
            B = Ah[:, indb]
            B_inv = invB(B_inv,p2,u)
           
            N = Ah[:, indn]
            iteration = iteration + 1
    return x[:n]
def MCF(A,b,L,P,K):
    S = A[2*L:,:P]
    R = A[:L,:P]
    indb = []
    [m,n] = A.shape
    u = b[:L,:];d = b[2*L:,:]
    
    A1 = np.hstack([S,np.zeros([K,L])])
    A2 = np.hstack([R,np.identity(L)])
    Ar = np.vstack([A1,A2])
    br = np.vstack([d,u])
    x = np.zeros([n,1])
    x[:P + L,:] = stage_one(Ar,br)
    geq = np.where(x[:P + L,:] > 0)[0]
    indb = np.append(indb,geq)
    #print(indb.size,K + L)
    if geq.size < K + L:
        radio = min(x[P:P + L,:]/u)
        M = 1.0 - radio# + 1e-2
        
        geq3 = np.arange(P + L,P + 2*L)
        indb = np.append(indb,geq3)
        if indb.size == K + 2*L - 1:
            pass
        else:
            er = K + 2*L - 1 - indb.size
            index = np.setdiff1d(np.arange(0,P + L),geq)
            ra = np.random.choice(index,er,replace = False)
            indb = np.append(indb,ra)
    else:
        radio = min(x[P:P + L,:]/u)
        ind = np.argmin(x[P:P + L,:]/u) + P + L
        geq3 = np.setdiff1d(np.arange(P + L,P + 2*L),ind)
        indb = np.append(indb,geq3)
        M = 1 - radio
    x[-1,0] = M
    indb = np.append(indb,P + 2*L)
    for i in range(L):
        x[P + L + i,0] = x[P + i,0] + (M - 1)*u[i,0]
    
    #print(indb.size)
    return x,indb
def stage_two(c,A,b,L,P,K):
    eps = 1e-6
    [m,n] = A.shape
    x,indb = MCF(A,b,L,P,K)
    #print(x.shape,A.shape,b.shape,indb.size,m)
    if np.linalg.norm(A@x - b) > eps or (any(x < 0) and min(x) < -eps):
        print('the init solution error')
        return 0
    else:
        print('the init solution come on')
    
    indb = indb.astype('int32')
    index = np.arange(0,n)
    indn = np.delete(index,indb)
    B = A[:,indb]
    
    B_inv = np.linalg.inv(B)
    
    N = A[:,indn]
    iteration = 0        # determine the max iteration
    
    while iteration < 3000:
        y = np.dot(B_inv.T, c[indb,:])
        s = c[indn,:] - N.T@y#C_N - N^TB^{-T}C_B
        value = (c*x).sum()
        if all(s >= 0) or (any(s < 0) and min(s) > -eps):#说明这个就是最优解
            print('the end epoch:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,np.linalg.norm(A@x - b),value))
            break
        z = np.random.rand(1)
        if (z >= 0 and z <= 1/3):
            q1, q2 = DanzigSelect(s, indn)
        elif (z >= 1/3 and z <= 2/3):
            A_re = B_inv@N
            A_n = np.linalg.norm(A_re,axis = 0,keepdims = True)
            q1, q2 = SteepSelect(A_n, s, indn)
        else:
            q1, q2 = Bland(s, indn)
        
        #q1是非基指标中需要被替换的索引
        u = B_inv@A[:, q1:(q1+1)]   
        if all(u <= 0):#说明最优值无界
            print('Linear program is unbounded.')
            break
            
        inds1 = np.where(u[:,0] > 0)[0]#找出所有大于0的索引
        inds2 = indb[inds1]#在基指标中把对应大于0的索引全找出来
        ratio = x[inds2,:]/u[inds1,:]
        x_q = min(ratio)        # a number
        p1 = inds2[ratio.argmin()]
        p2 = inds1[ratio.argmin()]
        
        x[indb,:] = x[indb,:]- x_q*u
        vec = np.zeros([indn.size,1])
        vec[q2,0] = x_q
        x[indn,:] = x[indn,:] + vec
        res = np.linalg.norm(A@x - b)
        if res > eps:
            print('blow up,epoch:%d,resid:%.2e'%(iteration,res))
            break
        if iteration%100 == 0:
            #print(np.linalg.det(B))
            print('epoch:%d,enter basis:%d,quit basis:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,q1,p1,res,value))
        indb = np.delete(indb,p2)
        indb = np.append(indb,q1)
            
        indn = np.delete(indn,q2)
        indn = np.append(indn,p1)
            
        B = A[:, indb]
        B_inv = invB(B_inv,p2,u)
           
        N = A[:, indn]
        iteration = iteration + 1
    return x




#L1 = 128
L1 = 64
P1 = 1024
k = 50
p = 0.02

np.random.seed(1)

# seed 1 with 512,1024, Danzig beats steepest
# seed 1 with 256, 512 steepest beats Danzig

R = np.random.binomial(1,p,[L1,P1])#二项分布，生成[L1,P1]行数组，元素全是0或1
d = np.random.gamma(2,2,[k,1])*1#gamma分布，
u = np.random.gamma(2,2,[L1,1])*3+20
s = P1 - k
ind = np.random.randint(0,k,s)#随机抽取0到K之间的整数
S = np.zeros((k,s))
for i in range(s):
    j = ind[i]
    S[j,i] = 1
S = np.hstack((np.identity(k),S))
n = P1+2*L1+1       #364
m = 2*L1+k          #196

# construct the vector c, only the last element is 1
c = np.zeros((n,1))
c[-1,0] = 1

# construct the matrix A
S1 = np.hstack((R,np.identity(L1),np.zeros((L1,L1)),np.zeros((L1,1))))
S2 = np.hstack((np.zeros((L1,P1)),np.identity(L1),-np.identity(L1),u))
S3 = np.hstack((S,np.zeros((k,(2*L1+1)))))
A = np.vstack((S1,S2,S3))


b = np.vstack((u,u,d))



#---------------------------------
L = R.shape[0];P = R.shape[1];K = S.shape[0]
m = 2*L + K;n = P + 2*L + 1
#print(type(R),type(S),type(u),type(d))
print(R.shape,S.shape,u.shape,d.shape)
print(max(R.reshape(-1,1)),min(R.reshape(-1,1)),max(S.reshape(-1,1)),min(S.reshape(-1,1)))
row1 = np.hstack([R,np.identity(L),np.zeros([L,L]),np.zeros([L,1])])
row2 = np.hstack([np.zeros([L,P]),np.identity(L),-np.identity(L),u])
row3 = np.hstack([S,np.zeros([K,2*L + 1])])
A = np.vstack([row1,row2,row3])
print(A.shape)
#print(np.linalg.matrix_rank(A),m)
b = np.vstack([u,u,d])
c = np.zeros([n,1]);c[-1,0] = 1.0
#piv = 'Bland'
#piv = 'Danzig'

st = time.time()
x = stage_two(c,A,b,L,P,K)
ela = time.time() - st

print(all(x >= 0),': use time:%.2f'%(ela))
geq = np.where(x > 0)[0]
eq = np.where(x == 0)[0]
leq = np.where(x < 0)[0]
print('geq num:%d,eq num:%d,leq num:%d, m=%d,n=%d'%(geq.size,eq.size,leq.size,m,n))
print('the nagetive element:')
print(x[leq])
#----------------------

singlesimplex.py

import numpy as np
import time
def Bland(s,indn):
    ind = np.where(s[:,0]<0)[0]
    q = ind[0]
    return indn[q], q

def DanzigSelect(s,indn):#indn非基索引
    q = s.argmin()
    return indn[q], q

def SteepSelect(A_n,s,indn):
    ratio = s/A_n.transpose()
    q = ratio.argmin()
    return indn[q], q
def invB(B_inv,p,u):
    Bv = B_inv.copy()
    Bv[p,:] = Bv[p,:]/u[p,0]
    m = u.size
    ind = []
    for i in range(m):
        if i == p:
            pass
        else:
            ind.append(i)
            Bv[i,:] = Bv[i,:] - u[i,0]*Bv[p,:]
    ind.append(p)
    Bv = Bv[ind,:]
    return Bv
def stage_one(A,b,piv):#使用大M方法求解初始基本可行解
    eps = 1e-6
    [m,n] = A.shape
    x = np.zeros([n + m,1]);x[n:,:] = b
    
    Ah = np.hstack([A,np.identity(m)])
    ch =  np.zeros([n + m,1]);ch[n:,:] = np.ones([m,1])
    indb = np.arange(m + n)[n:]#初始基指标索引
    index = np.arange(0,m + n)
    indn = np.delete(index,indb)
    B = Ah[:,indb]#初始B是单位矩阵
    B_inv = B
    
    N = Ah[:,indn]
    iteration = 0        # determine the max iteration
    
    while iteration < 1000:
        
        y = np.dot(B_inv.T, ch[indb,:])
        s = ch[indn,:] - N.T@y#C_N - N^TB^{-T}C_B
        if all(s >= 0):#说明这个就是最优解
            print('unbelievabel,the resid:%.2e, done,the epoch:%d'%(np.linalg.norm(Ah@x - b),iteration))
            break
        if piv == 'Danzig':
            q1, q2 = DanzigSelect(s, indn)
        if piv == 'steepest-edge':
            A_re = B_inv@N
            A_n = np.linalg.norm(A_re,axis = 0,keepdims = True)
            q1, q2 = SteepSelect(A_n, s, indn)
        if piv == 'Bland':
            q1, q2 = Bland(s, indn)
        #q1是非基指标中需要被替换的索引
        u = B_inv@Ah[:, q1:(q1+1)]   
        if all(u <= 0):#说明最优值无界
            print('Linear program is unbounded.')
            break
            
        inds1 = np.where(u[:,0] > 0)[0]#找出所有大于0的索引
        inds2 = indb[inds1]#在基指标中把对应大于0的索引全找出来
        ratio = x[inds2,:]/u[inds1,:]
        x_q = min(ratio)        # a number
        p1 = inds2[ratio.argmin()]
        p2 = inds1[ratio.argmin()]
        
        x[indb,:] = x[indb,:]- x_q*u
        vec = np.zeros([indn.size,1])
        vec[q2,0] = x_q
        x[indn,:] = x[indn,:] + vec
        aug = np.linalg.norm(Ah@x - b)
        if aug > eps:
            print('blow up,epoch:%d,aug resid:%.2e'%(iteration,aug))
            break
        value = x[n:,:].sum()
        if iteration%100 == 0:
            
            print('epoch:%d,enter basis:%d,quit basis:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,q1,p1,aug,value))
        
        if value < eps and all(s) >= 0:
            print('the end epoch:%d,opt val:%.2e,the real resid:%.2e'
                  %(iteration,value,np.linalg.norm(A@x[:n,:] - b)))
            break
        else:
            
            indb = np.delete(indb,p2)
            indb = np.append(indb,q1)
            
            indn = np.delete(indn,q2)
            indn = np.append(indn,p1)
            
            B = Ah[:, indb]
            B_inv = invB(B_inv,p2,u)
           
            N = Ah[:, indn]
            iteration = iteration + 1
    return x[:n]
def MCF(A,b,L,P,K):
    S = A[2*L:,:P]
    R = A[:L,:P]
    indb = []
    [m,n] = A.shape
    u = b[:L,:];d = b[2*L:,:]
    
    A1 = np.hstack([S,np.zeros([K,L])])
    A2 = np.hstack([R,np.identity(L)])
    Ar = np.vstack([A1,A2])
    br = np.vstack([d,u])
    x = np.zeros([n,1])
    x[:P + L,:] = stage_one(Ar,br,piv = 'Bland')
    geq = np.where(x[:P + L,:] > 0)[0]
    indb = np.append(indb,geq)
    #print(indb.size,K + L)
    if geq.size < K + L:
        radio = min(x[P:P + L,:]/u)
        M = 1.0 - radio# + 1e-2
        
        geq3 = np.arange(P + L,P + 2*L)
        indb = np.append(indb,geq3)
        if indb.size == K + 2*L - 1:
            pass
        else:
            er = K + 2*L - 1 - indb.size
            index = np.setdiff1d(np.arange(0,P + L),geq)
            ra = np.random.choice(index,er,replace = False)
            indb = np.append(indb,ra)
    else:
        radio = min(x[P:P + L,:]/u)
        ind = np.argmin(x[P:P + L,:]/u) + P + L
        geq3 = np.setdiff1d(np.arange(P + L,P + 2*L),ind)
        indb = np.append(indb,geq3)
        M = 1 - radio
    x[-1,0] = M
    indb = np.append(indb,P + 2*L)
    for i in range(L):
        x[P + L + i,0] = x[P + i,0] + (M - 1)*u[i,0]
    
    #print(indb.size)
    return x,indb
def stage_two(c,A,b,piv,L,P,K):
    eps = 1e-6
    [m,n] = A.shape
    x,indb = MCF(A,b,L,P,K)
    #print(x.shape,A.shape,b.shape,indb.size,m)
    if np.linalg.norm(A@x - b) > eps or (any(x < 0) and min(x) < -eps):
        print('the init solution error')
        return 0
    else:
        print('the init solution come on')
    
    indb = indb.astype('int32')
    index = np.arange(0,n)
    indn = np.delete(index,indb)
    B = A[:,indb]
    
    B_inv = np.linalg.inv(B)
    
    N = A[:,indn]
    iteration = 0        # determine the max iteration
    
    while iteration < 2000:
        y = np.dot(B_inv.T, c[indb,:])
        s = c[indn,:] - N.T@y#C_N - N^TB^{-T}C_B
        value = (c*x).sum()
        if all(s >= 0) or (any(s < 0) and min(s) > -eps):#说明这个就是最优解
            print('the end epoch:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,np.linalg.norm(A@x - b),value))
            break
        if piv == 'Danzig':
            q1, q2 = DanzigSelect(s, indn)
        if piv == 'steepest-edge':
            A_re = B_inv@N
            A_n = np.linalg.norm(A_re,axis = 0,keepdims = True)
            q1, q2 = SteepSelect(A_n, s, indn)
        if piv == 'Bland':
            q1, q2 = Bland(s, indn)
        #q1是非基指标中需要被替换的索引
        u = B_inv@A[:, q1:(q1+1)]   
        if all(u <= 0):#说明最优值无界
            print('Linear program is unbounded.')
            break
            
        inds1 = np.where(u[:,0] > 0)[0]#找出所有大于0的索引
        inds2 = indb[inds1]#在基指标中把对应大于0的索引全找出来
        ratio = x[inds2,:]/u[inds1,:]
        x_q = min(ratio)        # a number
        p1 = inds2[ratio.argmin()]
        p2 = inds1[ratio.argmin()]
        
        x[indb,:] = x[indb,:]- x_q*u
        vec = np.zeros([indn.size,1])
        vec[q2,0] = x_q
        x[indn,:] = x[indn,:] + vec
        res = np.linalg.norm(A@x - b)
        if res > eps:
            print('blow up,epoch:%d,resid:%.2e'%(iteration,res))
            break
        if iteration%100 == 0:
            #print(np.linalg.det(B))
            print('epoch:%d,enter basis:%d,quit basis:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,q1,p1,res,value))
        indb = np.delete(indb,p2)
        indb = np.append(indb,q1)
            
        indn = np.delete(indn,q2)
        indn = np.append(indn,p1)
            
        B = A[:, indb]
        B_inv = invB(B_inv,p2,u)
           
        N = A[:, indn]
        iteration = iteration + 1
    return x




L1 = 128
P1 = 1024
k = 50
p = 0.02

np.random.seed(1)

# seed 1 with 512,1024, Danzig beats steepest
# seed 1 with 256, 512 steepest beats Danzig

R = np.random.binomial(1,p,[L1,P1])#二项分布，生成[L1,P1]行数组，元素全是0或1
d = np.random.gamma(2,2,[k,1])*1#gamma分布，
u = np.random.gamma(2,2,[L1,1])*3+20
s = P1 - k
ind = np.random.randint(0,k,s)#随机抽取0到K之间的整数
S = np.zeros((k,s))
for i in range(s):
    j = ind[i]
    S[j,i] = 1
S = np.hstack((np.identity(k),S))
n = P1+2*L1+1       #364
m = 2*L1+k          #196

# construct the vector c, only the last element is 1
c = np.zeros((n,1))
c[-1,0] = 1

# construct the matrix A
S1 = np.hstack((R,np.identity(L1),np.zeros((L1,L1)),np.zeros((L1,1))))
S2 = np.hstack((np.zeros((L1,P1)),np.identity(L1),-np.identity(L1),u))
S3 = np.hstack((S,np.zeros((k,(2*L1+1)))))
A = np.vstack((S1,S2,S3))


b = np.vstack((u,u,d))



#---------------------------------
L = R.shape[0];P = R.shape[1];K = S.shape[0]
m = 2*L + K;n = P + 2*L + 1
#print(type(R),type(S),type(u),type(d))
print(R.shape,S.shape,u.shape,d.shape)
print(max(R.reshape(-1,1)),min(R.reshape(-1,1)),max(S.reshape(-1,1)),min(S.reshape(-1,1)))
row1 = np.hstack([R,np.identity(L),np.zeros([L,L]),np.zeros([L,1])])
row2 = np.hstack([np.zeros([L,P]),np.identity(L),-np.identity(L),u])
row3 = np.hstack([S,np.zeros([K,2*L + 1])])
A = np.vstack([row1,row2,row3])
#print(np.linalg.matrix_rank(A),m)
b = np.vstack([u,u,d])
c = np.zeros([n,1]);c[-1,0] = 1.0
#piv = 'Bland'
piv = 'Danzig'

st = time.time()
x = stage_two(c,A,b,piv,L,P,K)
ela = time.time() - st

print(all(x >= 0),': use time:%.2f'%(ela))
geq = np.where(x > 0)[0]
eq = np.where(x == 0)[0]
leq = np.where(x < 0)[0]
print('geq num:%d,eq num:%d,leq num:%d, m=%d,n=%d'%(geq.size,eq.size,leq.size,m,n))
print('the nagetive element:')
print(x[leq])
#----------------------

在这里插入图片描述

multisimplex.py

在这里插入图片描述

import numpy as np
import time
def Bland(s,indn):
    ind = np.where(s[:,0]<0)[0]
    q = ind[0]
    return indn[q], q

def DanzigSelect(s,indn):#indn非基索引
    q = s.argmin()
    return indn[q], q

def SteepSelect(A_n,s,indn):
    ratio = s/A_n.transpose()
    q = ratio.argmin()
    return indn[q], q
def invB(B_inv,p,u):
    Bv = B_inv.copy()
    Bv[p,:] = Bv[p,:]/u[p,0]
    m = u.size
    ind = []
    for i in range(m):
        if i == p:
            pass
        else:
            ind.append(i)
            Bv[i,:] = Bv[i,:] - u[i,0]*Bv[p,:]
    ind.append(p)
    Bv = Bv[ind,:]
    return Bv
def stage_one(A,b,piv):#使用大M方法求解初始基本可行解
    eps = 1e-6
    [m,n] = A.shape
    x = np.zeros([n + m,1]);x[n:,:] = b
    
    Ah = np.hstack([A,np.identity(m)])
    ch =  np.zeros([n + m,1]);ch[n:,:] = np.ones([m,1])
    indb = np.arange(m + n)[n:]#初始基指标索引
    index = np.arange(0,m + n)
    indn = np.delete(index,indb)
    B = Ah[:,indb]#初始B是单位矩阵
    B_inv = B
    
    N = Ah[:,indn]
    iteration = 0        # determine the max iteration
    
    while iteration < 1000:
        
        y = np.dot(B_inv.T, ch[indb,:])
        s = ch[indn,:] - N.T@y#C_N - N^TB^{-T}C_B
        if all(s >= 0):#说明这个就是最优解
            print('unbelievabel,the resid:%.2e, done,the epoch:%d'%(np.linalg.norm(Ah@x - b),iteration))
            break
        if piv == 'Danzig':
            q1, q2 = DanzigSelect(s, indn)
        if piv == 'steepest-edge':
            A_re = B_inv@N
            A_n = np.linalg.norm(A_re,axis = 0,keepdims = True)
            q1, q2 = SteepSelect(A_n, s, indn)
        if piv == 'Bland':
            q1, q2 = Bland(s, indn)
        #q1是非基指标中需要被替换的索引
        u = B_inv@Ah[:, q1:(q1+1)]   
        if all(u <= 0):#说明最优值无界
            print('Linear program is unbounded.')
            break
            
        inds1 = np.where(u[:,0] > 0)[0]#找出所有大于0的索引
        inds2 = indb[inds1]#在基指标中把对应大于0的索引全找出来
        ratio = x[inds2,:]/u[inds1,:]
        x_q = min(ratio)        # a number
        p1 = inds2[ratio.argmin()]
        p2 = inds1[ratio.argmin()]
        
        x[indb,:] = x[indb,:]- x_q*u
        vec = np.zeros([indn.size,1])
        vec[q2,0] = x_q
        x[indn,:] = x[indn,:] + vec
        aug = np.linalg.norm(Ah@x - b)
        if aug > eps:
            print('blow up,epoch:%d,aug resid:%.2e'%(iteration,aug))
            break
        value = x[n:,:].sum()
        if iteration%100 == 0:
            
            print('epoch:%d,enter basis:%d,quit basis:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,q1,p1,aug,value))
        
        if value < eps:
            print('the end epoch:%d,opt val:%.2e,the real resid:%.2e'
                  %(iteration,value,np.linalg.norm(A@x[:n,:] - b)))
            break
        else:
            
            indb = np.delete(indb,p2)
            indb = np.append(indb,q1)
            
            indn = np.delete(indn,q2)
            indn = np.append(indn,p1)
            
            B = Ah[:, indb]
            B_inv = invB(B_inv,p2,u)
           
            N = Ah[:, indn]
            iteration = iteration + 1
    return x[:n,:]
def MCF(A,b,L,P,K):
    S = A[2*L:,:P]
    R = A[:L,:P]
    indb = []
    [m,n] = A.shape
    u = b[:L,:];d = b[2*L:,:]
    
    A1 = np.hstack([S,np.zeros([K,L])])
    A2 = np.hstack([R,np.identity(L)])
    Ar = np.vstack([A1,A2])
    br = np.vstack([d,u])
    x = np.zeros([n,1])
    x[:P + L,:] = stage_one(Ar,br,piv = 'Bland')
    geq = np.where(x[:P + L,:] > 0)[0]
    indb = np.append(indb,geq)
    #print(indb.size,K + L)
    if geq.size < K + L:
        radio = min(x[P:P + L,:]/u)
        M = 1.0 - radio# + 1e-2
        
        geq3 = np.arange(P + L,P + 2*L)
        indb = np.append(indb,geq3)
        if indb.size == K + 2*L - 1:
            pass
        else:
            er = K + 2*L - 1 - indb.size
            index = np.setdiff1d(np.arange(0,P + L),geq)
            ra = np.random.choice(index,er,replace = False)
            indb = np.append(indb,ra)
    else:
        radio = min(x[P:P + L,:]/u)
        ind = np.argmin(x[P:P + L,:]/u) + P + L
        geq3 = np.setdiff1d(np.arange(P + L,P + 2*L),ind)
        indb = np.append(indb,geq3)
        M = 1 - radio
    x[-1,0] = M
    indb = np.append(indb,P + 2*L)
    for i in range(L):
        x[P + L + i,0] = x[P + i,0] + (M - 1)*u[i,0]
    
    #print(indb.size)
    return x,indb
def stage_two(c,A,b,L,P,K):
    eps = 1e-6
    [m,n] = A.shape
    x,indb = MCF(A,b,L,P,K)
    
    if np.linalg.norm(A@x - b) > eps or (any(x < 0) and min(x) < -eps):
        print('the init solution error')
        return 0
    else:
        print('the init solution come on')
    
    indb = indb.astype('int32')
    index = np.arange(0,n)
    indn = np.delete(index,indb)
    B = A[:,indb]
    
    B_inv = np.linalg.inv(B)
    
    N = A[:,indn]
    iteration = 0        # determine the max iteration
    value = 0
    while iteration < 2000:
        y = np.dot(B_inv.T, c[indb,:])
        s = c[indn,:] - N.T@y#C_N - N^TB^{-T}C_B
        
        if all(s >= 0) or (any(s < 0) and min(s) > -eps):#说明这个就是最优解
            print('the end epoch:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,np.linalg.norm(A@x - b),value))
            break
        xd = x.copy()
        xs = x.copy()
        xb = x.copy()
        allx = []
        allq = []
        allp = []
        allu = []
        allval = []
        #------------------------
        q1, q2 = DanzigSelect(s, indn)
        allq.append(q1)
        allq.append(q2)
        u = B_inv@A[:, q1:(q1+1)]
        inds1 = np.where(u[:,0] > 0)[0]#找出所有大于0的索引
        inds2 = indb[inds1]#在基指标中把对应大于0的索引全找出来
        ratio = xd[inds2,:]/u[inds1,:]
        x_q = min(ratio)        # a number
        p1 = inds2[ratio.argmin()]
        p2 = inds1[ratio.argmin()]
        allp.append(p1)
        allp.append(p2)
        xd[indb,:] = xd[indb,:]- x_q*u
        vec = np.zeros([indn.size,1])
        vec[q2,0] = x_q
        xd[indn,:] = xd[indn,:] + vec
        
        val = (c*xd).sum()
        
        allx.append(xd)
        
        allu.append(u)
        allval.append(val)
        #------------------------
        A_re = B_inv@N
        A_n = np.linalg.norm(A_re,axis = 0,keepdims = True)
        q1, q2 = SteepSelect(A_n, s, indn)
        allq.append(q1)
        allq.append(q2)
        u = B_inv@A[:, q1:(q1+1)]
        inds1 = np.where(u[:,0] > 0)[0]#找出所有大于0的索引
        inds2 = indb[inds1]#在基指标中把对应大于0的索引全找出来
        ratio = xs[inds2,:]/u[inds1,:]
        x_q = min(ratio)        # a number
        p1 = inds2[ratio.argmin()]
        p2 = inds1[ratio.argmin()]
        allp.append(p1)
        allp.append(p2)
        xs[indb,:] = xs[indb,:]- x_q*u
        vec = np.zeros([indn.size,1])
        vec[q2,0] = x_q
        xs[indn,:] = xs[indn,:] + vec
        
        val = (c*xs).sum()
        
        allx.append(xs)
        
        allu.append(u)
        allval.append(val)
        #-------------------------
        q1, q2 = Bland(s, indn)
        allq.append(q1)
        allq.append(q2)
        u = B_inv@A[:, q1:(q1+1)]
        inds1 = np.where(u[:,0] > 0)[0]#找出所有大于0的索引
        inds2 = indb[inds1]#在基指标中把对应大于0的索引全找出来
        ratio = xb[inds2,:]/u[inds1,:]
        x_q = min(ratio)        # a number
        p1 = inds2[ratio.argmin()]
        p2 = inds1[ratio.argmin()]
        allp.append(p1)
        allp.append(p2)
        xb[indb,:] = xb[indb,:]- x_q*u
        vec = np.zeros([indn.size,1])
        vec[q2,0] = x_q
        xb[indn,:] = xb[indn,:] + vec
        
        val = (c*xb).sum()
        
        allx.append(xb)
        
        allu.append(u)
        allval.append(val)
        #----------------
        mind = np.argmin(allval)
        u = allu[mind]
        x = allx[mind]
        
        value = allval[mind]
        q1,q2 = allq[2*mind],allq[2*mind + 1]
        p1,p2 = allp[2*mind],allp[2*mind + 1]
        allu.clear();allx.clear();allval.clear();allq.clear();allp.clear()
        #q1是非基指标中需要被替换的索引
        if all(u <= 0):#说明最优值无界
            print('Linear program is unbounded.')
            break
        res = np.linalg.norm(A@x - b)
        if res > eps:
            print('blow up,epoch:%d,resid:%.2e'%(iteration,res))
            break
        if iteration%100 == 0:
            #print(np.linalg.det(B))
            print('epoch:%d,enter basis:%d,quit basis:%d,resid:%.2e,opt val:%.2e'
                  %(iteration,q1,p1,res,value))
        indb = np.delete(indb,p2)
        indb = np.append(indb,q1)
            
        indn = np.delete(indn,q2)
        indn = np.append(indn,p1)
            
        B = A[:, indb]
        B_inv = invB(B_inv,p2,u)
           
        N = A[:, indn]
        iteration = iteration + 1
    return x






L1 = 128
P1 = 1024
k = 50
p = 0.02

np.random.seed(1)

# seed 1 with 512,1024, Danzig beats steepest
# seed 1 with 256, 512 steepest beats Danzig

R = np.random.binomial(1,p,[L1,P1])#二项分布，生成[L1,P1]行数组，元素全是0或1
d = np.random.gamma(2,2,[k,1])*1#gamma分布，
u = np.random.gamma(2,2,[L1,1])*3+20
s = P1 - k
ind = np.random.randint(0,k,s)#随机抽取0到K之间的整数
S = np.zeros((k,s))
for i in range(s):
    j = ind[i]
    S[j,i] = 1
S = np.hstack((np.identity(k),S))
n = P1+2*L1+1       #364
m = 2*L1+k          #196

# construct the vector c, only the last element is 1
c = np.zeros((n,1))
c[-1,0] = 1

# construct the matrix A
S1 = np.hstack((R,np.identity(L1),np.zeros((L1,L1)),np.zeros((L1,1))))
S2 = np.hstack((np.zeros((L1,P1)),np.identity(L1),-np.identity(L1),u))
S3 = np.hstack((S,np.zeros((k,(2*L1+1)))))
A = np.vstack((S1,S2,S3))


b = np.vstack((u,u,d))


#---------------------------------
L = R.shape[0];P = R.shape[1];K = S.shape[0]
m = 2*L + K;n = P + 2*L + 1
#print(type(R),type(S),type(u),type(d))
print(R.shape,S.shape,u.shape,d.shape)
print(max(R.reshape(-1,1)),min(R.reshape(-1,1)),max(S.reshape(-1,1)),min(S.reshape(-1,1)))
row1 = np.hstack([R,np.identity(L),np.zeros([L,L]),np.zeros([L,1])])
row2 = np.hstack([np.zeros([L,P]),np.identity(L),-np.identity(L),u])
row3 = np.hstack([S,np.zeros([K,2*L + 1])])
A = np.vstack([row1,row2,row3])
#print(np.linalg.matrix_rank(A),m)
b = np.vstack([u,u,d])
c = np.zeros([n,1]);c[-1,0] = 1.0

st = time.time()
x = stage_two(c,A,b,L,P,K)
ela = time.time() - st

print(all(x >= 0),': use time:%.2f'%(ela))
geq = np.where(x > 0)[0]
eq = np.where(x == 0)[0]
leq = np.where(x < 0)[0]
print('geq num:%d,eq num:%d,leq num:%d, m=%d,n=%d'%(geq.size,eq.size,leq.size,m,n))
#print(geq.size,eq.size,leq.size,m)
print('the nagetive element:')
print(x[leq])
#----------------------

Galerkin码农选手

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
打赏
0
评论
单纯形法求解线性规划

问题描述考虑从源节点到端节点传输数据的网络，我们首先引入一些定义：\∙\bullet∙ K\quad flows：传输K中不同的数据。\∙\bullet∙ L\quad link：连接网络节点的L条无向边。\∙\bullet∙ clc_lcl：在第lll条边上传输数据的最大容量。\∙\bullet∙ PkP_kPk：对于不同的k\quad flow，运输成功可选择的路path的数目。\∙\bullet∙ xk=(xk,1,…,xk,PK)T\mathbf{x_k}=\left(x_{k,1}
复制链接

扫一扫