python 线性规划求解

参考文献：《基于线性规划方法处理c波段双线偏振 多普勒天气雷达差分相位"c数据》

最后求解不出来，请大佬指点



########
#加载所需包
import math
import time
# from itkwidgets import view
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np
from matplotlib import colors
from scipy.interpolate import griddata
import scipy.misc
import cartopy.crs as ccrs
import cartopy.feature as cfeat
import cinrad
import xarray as xr
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
from cinrad.io import CinradReader, StandardData
from cinrad.visualize import PPI, Section
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd
import plotly.express as px
from metpy.io import Level3File
from metpy.plots import add_metpy_logo,colortables
%matplotlib inline
# 注意:x是一维序列,r是一个数字
#import pulp as pp #####线性规划
#from pulp import LpVariable, LpProblem, lpSum, LpMinimize, LpStatus, value
import  pulp as pp
####################################################################################################   

def printMatrix(matrix):   #####显示矩阵
    for row in  matrix:
        print( ' '.join( str(i) for i in row))
###############################################读取数据###################################################      
f = cinrad.io.StandardData(r"G:\2022\Z9092.20220509.000135.AR2")
# 查看包含产品类型
f.available_product(0)
#雷达扫描仰角
f.available_tilt("PHI")  ###差分相移
#f.available_tilt("PHI")  ###差分相移
ele = 0 #选择第1个仰角
radius = 100 #绘制图像的范围大小，单位km
r = f.get_data(ele,radius,"PHI") #选择差分相移数据
rho = f.get_data(ele,radius,"RHO") ##相关系数 只有100km
#print(r.dims)
#print(r.variables)
#print(r.coords)
az = np.round(np.rad2deg(r.azimuth)) #####  将角度从弧度转换为度数后 四舍五入    az方位角
#print(az)
r['azimuth']=az ###改变r的坐标轴数值  az方位角
rho['azimuth']=az

ref=r.PHI     ######提取需要的变量
rrho=rho.RHO   ###azimuth: 361distance: 400
################################线性规划  #########################################################   
nshu=len(ref[0,0:20])  #### bi有n个数据
m=nshu-4    #####5点差分滤波

bi=np.array(ref[0,0:20])  ##转化为矩阵
print("差分相位bi为",bi)
print(type(bi))  ###确
# 目标函数的系数
c = [1 for i in range(nshu)]+[0 for i in range(nshu)]
print("目标函数的系数c",c)
print(type(c))   ###确
#C=np.array(c)

Mn=[None]*m
Zn4=np.array([0]*m)
print("中间变量Zn4",Zn4)
print(type(Zn4))  ###


for i in range(m):    ##创建一个n*m的零矩阵
    Mn=[0]*nshu
i=0
Zn=np.array(Mn)   ###Zn 为零矩阵 （n-4,n）
for j in range(nshu):
   # print(j)
    if i<=m-1:
        Mn[j]=-0.2
        Mn[j][i+1]=-0.1
        Mn[j][i+2]=0.0
        Mn[j][i+3]=0.1
        Mn[j][i+4]=0.2
        i=i+1
print("变量Mn",Mn)  
M=np.array(Mn)
print("array变量M",M[1][1])###微分滤波器Mn （n-4,n）
print(type(M))
print("array变量M",M)
printMatrix(M)
print(M.shape)### 可以得到数组（或矩阵）的行数和列数(m,n)
print(Zn.shape)
print("零矩阵Zn",Zn)
print(type(Zn))


I=np.eye(nshu)
print("单位变量I",I)
print(type(I))


A1=np.append(I,-I,axis=1)
print("对应变量A1",A1)
#printMatrix(A1)
print(type(A1))
print(A1.shape)
A2=np.append(I,I,axis=1)
A3=np.append(Zn,M,axis=1)
print("对应变量A3",A3)
#printMatrix(A1)
print(type(A3))
print(A3.shape)
A4=np.append(A1,A2,axis=0)
print("对应变量A4",A4)
#printMatrix(A1)
print(type(A4))
print(A4.shape)

A = np.append(A4,A3,axis=0)######np.array([[I, -I],[I,I],[Zn,M]],dtype=object)   ######需
print("对应变量A",A)
print(type(A))
print(A.shape)

B1= np.append(-bi,bi)####np.array([-bi,bi,Zn4],dtype=object)
B= np.append(B1,Zn4)
print("对应变量B",B)
print(type(B))
print(B.shape)


###################################################################################################
# 确定最大最小化问题，当前确定的是最小化问题  
#m = pp.LpProblem(sense=pp.LpMinimize)
prob = pp.LpProblem('Ensemble', LpMinimize)

# 定义n个变量放到列表中
x = [pp.LpVariable(f'x{i}', cat = 'Continuous') for i in range(nshu)]
print("求解变量",x)
print(type(x))  ##

# 定义n个变量放到列表中
z = [pp.LpVariable(f'z{i}', cat = 'Continuous') for i in range(nshu)]
print("中间变量",z)
print(type(z))  ##
print(z[0])


xc=np.append(z,x)##np.array([z,x],dtype=object)
print("对应变量xc",xc)
print(type(xc))
print(xc.shape)

prob += ( pp.lpDot(c, xc))  # 设置目标函数 f(x)
prob += ( pp.lpDot(A, xc)>= B)  # 不等式约束

# 设置比较条件
for i in range(nshu):
    prob += ((z-x) >= -bi)

# 设置比较条件
for i in range(nshu):
    prob += ((z+x) >= bi)


prob.solve()
print(prob.name)  # 输出求解状态
# 输出结果
print(f'优化结果：{pp.value(prob.objective)}')
print(f'参数取值：{[pp.value(var) for var in x]}')