【数学建模基础】约束性拟合（线性）--python_cp.variable()函数功能

前言：        

        线性约束性拟合分为拟合与约束两部分，即对‘拟合’的估计参数进行‘约束’（本文运用最小二乘法进行解题。）例如:

已知x,y的观测值如下表所示。用给定数据拟合函数,且满足

 分析：

         可知为待估参数，为约束条件。

import numpy as np
import cvxpy as cp
x=np.array([3,5,6,7,4,8,5,9])#必须将列表转化为数组才能进行下一步计算
y=np.array([4,9,5,3,8,5,8,5])
#含义为预拟合2个参数,返回到一个t数组中并且这两个参数的数值都要大于0(pos=True)
t=cp.Variable(2,pos=True)
#对所给列表中数组先按列进行拼接形成一个2*5的矩阵然后进行转置|a*exp(x)+b*ln(x)
c=np.vstack([np.exp(x),np.log(x)]).T
#min sum(a*exp(x)+b*ln(x)-y)^2
obj=cp.Minimize(cp.sum_squares(c@ t -y0))
#约束条件
con=[sum(t)<=1]
prob=cp.Problem(obj,con)
prob.solve(solver='CVXOPT')
print('1',t.value)

 结果为：[4.11416418e-04 9.99588813e-01]（拟合数）

 异同：

        上述代码与规划密切相关，属于‘线性规划’方法，与‘拟合’方法毫无关系，但是可以通过此来方法解决线性约束拟合问题，下述规划问题代码与上述代码极其相似。

例: 为了生产需要，某工厂一条生产线需要24h运转，但是每天不同时间段所需要的工人最低数量不相同，具体如下表所示。已知每名工人连续工作时间为8h。则该工厂应该为生产线配置多少名工人才能保证生产线正常运转 

import numpy as np
import cvxpy as cp
a=np.array([35,40,50,45,55,30])
x=cp.Variable(6,integer=True)#作用于待估参数
obj1=cp.Minimize(sum(x))
cons1=[x>=0]
for i in range(6):
    cons1.append(x[(i-1)%6]+x[i]>=a[i])
prob1=cp.Problem(obj1,cons1)
prob1.solve(solver='GLPK_MI')
print(prob1.value,x.value,sep='\n')

 结果：140(目标最小数总数).[35,5,45,25,30, 0]（格未知数最佳值）

线性最小二乘法 '拟合'常规方法如下：

        例：某天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，两坐标轴上的单位长度取为1天文测量单位。在5个不同时间对小行星做了5次观测，测得轨道上5个点的坐标数据如下表所示。由开普勒第一定律可知，小行星的轨道为一椭圆，一般方程为：

$a_{\1}x^{\2}+a_{\2}xy+a_{\3}y^{\2}+a_{\4}x+a_{\5}y+1=0.$

据此建立行星运行轨道方程，并画出轨道曲线。

import numpy as np
import matplotlib.pyplot as plt
x=np.array([5.764,6.286,6.759,7.168,7.408])#必须是数组才能进行计算
y=np.array([0.648,1.202,1.823,2.526,3.360])
A=np.vstack([x**2,x*y,y**2,x,y]).T#等式两边给不能是零！
z=-np.ones(len(A))
p=np.linalg.pinv(A) @ z
print(np.round(p,4))
x=np.linspace(3,8,100)
y=np.linspace(-1,5,100)
x,y=np.meshgrid(x,y)
f=lambda x,y:p[0]*x**2+p[1]*x*y+p[2]*y**2+p[3]*x+p[4]*y
z=f(x,y)
plt.contour(x,y,z,[-1])
plt.show()

结果为：$a_{\1}=0.0508,a_{\2}=-0.0702,a_{\3}=0.0381,a_{\4}=0.2643.$

 

 方法中的关键函数：

        规划： 

                cvxpy.Variable():可以指定拟合参数的个数，并且可以对参数进行控制,例如是否为整数，是否大于零等等。

               cvxpy.Minimize（_@_）和cvxpy.Maximize(_@_):将需要规划的未知量与其参数结合形成目标函数object，其中cvxpy.sum_square(_@_)等同于$\sum \left |x_{\i} \right |^{\2}$

                 cvxpy.Problem(object,constraints):单目标规划，将object（目标函数）与constraints（约束条件）对应并求解

        拟合:

                numpy.vstack()和numpy.hstack():给定的参数必须是数组或者是元组，并且数组和元组中的元素的维数相同，对于元素都是一维的数组来说，numpy.vstack作用为将各元素对应形成一个矩阵，大小为（给定列表或元组的长度，各元素的维度），对于元素都是一维的数组来说，numpy.hstack作用为将各堆叠形成一个数组，大小为（1，各元素维数之和）。

                注：列表中的各元素中的元素个数和列表中元素的维度要求一致

import numpy as np
a=np.array([1,2,3])
b=np.array([7,8,9])
c=np.vstack([a,b])
d=np.hstack([a,b])
print(a,b,c,d,sep='\n')

         numpy.linalg.pinv() @_:将所给参数进行拟合，可以视作一个等式，而且等式两边都不能为零否则拟合结果会得到异常值。此函数也可以由numpy.polyfit和numpy.polyval(计算函数值)代替

 来源：

《python数学建模算法于应用》----司守奎，孙玺菁。