非线性规划、例题的matlab\python实现、选址问题的matlab实现

非线性规划模型

• 目标函数或约束条件中包含非线性函数的数学规划问题称为非线性规划问题

$$\begin{gathered} \min f(\mathbf{x}) \\ s.t.\{\begin{array}{l}A\cdot \mathbf{x}\le \mathbf{b}\\ Aeq\cdot \mathbf{x}=beq\\ c(\mathbf{x})\le 0\\ ceq(\mathbf{x})=0\\ lb\le \mathbf{x}\le ub\end{array} \end{gathered}$$

• Matlab求解非线性规划问题的函数

[x,fval]=fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlfun,option)
fun:目标函数的M文件名
nonlfun:非线性约束条件函数的M文件名
option:可以给定求解的算法，可以使用默认参数设置

• 注意

1. 非线性规划问题对于初始值x0的选取很重要，因为费霞凝规划的算法求解出来的是一个局部最优解
2. “全局最优解”的两种思路：

①给定不同的初值，在里面找到一个最优解
②先用蒙特卡洛模拟得到一个解，然后将这个解作为初始值来求最优解（推荐）

3. option选项有四种算法：interior-point(内点法)、sqp(序列二次规划法)、action-set(有效集法)、trust-region-reflective(信赖域反射算法)，不同的算法有各自的优缺点，我们可以该边求解的算法来看求解的结果是否更优

例题并求解

$$\begin{gathered} \min f(\mathbf{x})=x_1^2+x_2^2+x_3^2+8 \\ s.t.\{\begin{array}{l}{x}_1^2-x_2+x_3^2\ge 0\\ x_1+x_2^2+x_3^3\le 20\\ -x_1-x_2^2+2=0\\ x_2+2x_3^2=3\\ x_1,x_2,x_3\ge 0\end{array} \end{gathered}$$

• Matlab求解

%% fun2.m定义目标函数
function f = fun2(x)
f=sum(x.^2)+8;
end

%% nonlfun2.m定义非线性约束
function [g,h]=nonlfun2(x)
%g 非线性不等式约束；h 非线性等式约束
g=[-x(1)^2+x(2)-x(3)^2;x(1)+x(2)^2+x(3)^3]-20;
h=[-x(1)-x(2)^2+2;x(2)+2*x(3)^2-3];
end

%% 主程序
clc;clear;
%% 用蒙特卡罗方法找初始值
f_min=+inf;
x0=[0 ;0 ;0];
x1=unifrnd(0,2,10^6,1);
x2=unifrnd(0,1.5,10^6,1);
x3=unifrnd(0,1.3,10^6,1);
tic
for i=1:10^6
    x=[x1(i),x2(i),x3(i)];
    f=fun2(x);
    [g,h]=nonlfun2(x);
    if isempty(find(g>0)) & abs(h(1))<0.01 & abs(h(2))<0.01
        if f<f_min
            x0=x;f_min=f;
        end
    end
end
x0
f_min
toc
%% 代入求得的初值求解
[x,fval]=fmincon(@fun2,[0.3471,1.2852,0.9278],[],[],[],[],zeros(3,1),[],@nonlfun2)

结果为：x = 0.4011 1.2645 0.9315 fval = 10.6275
对比《数学建模算法与应用》上的结果，此结果更优

• python求解

#目标函数
fun=lambda x:x[0]**2+x[1]**2+x[2]**2+8
# 'eq'表示表达式等于0，'ineq'表示表达式大于等于0
constraints=[{'type': 'ineq', 'fun': lambda x:x[0]**2 - x[1] + x[2]**2},
             {'type': 'ineq', 'fun': lambda x: 20 - x[0] - x[1]**2 - x[2]**3},
             {'type': 'eq', 'fun': lambda x: -x[0] - x[1]**2 +2},
             {'type': 'eq', 'fun': lambda x: x[1] + 2*x[2]**2 - 3}]
#上下界
bounds = [(0,None)]*3

#蒙特卡洛法设置初始值
f_min=float('inf')
x0=[0 ,0 ,0]
x1=2*np.random.rand(10**6,1)
x2=1.5*np.random.rand(10**6,1)
x3=1.3*np.random.rand(10**6,1)
for i in range(10**6):
    x=[x1[i],x2[i],x3[i]]
    x=[x[0][0],x[1][0],x[2][0]]
    f=fun(x)
    g=np.array([constraints[0]['fun'](x),constraints[1]['fun'](x)])
    h=np.array([constraints[2]['fun'](x),constraints[3]['fun'](x)])
    if all(g)>=0 & all(abs(h))<0.00001:
        x0=x
        f_min=f
print(x0)
print(f_min)

#求解
res = minimize(fun, x0, method='SLSQP',constraints=constraints)
res

结果：fun: 10.651091840570237
x: array([0.55216734, 1.20325918, 0.94782404])

非线性规划的典型例题——选址问题

• 题目：临时料场A（5，1），B（2，7），日储量各20吨。

工地位置坐标及日需求量
	1	2	3	4	5	6
横坐标	1.25	8.75	0.5	5.75	3	7.25
纵坐标	1.25	0.75	4.75	5	6.5	7.25
日需求量	3	5	4	7	6	11

（1）试制定每天的供应计划，即从两料场分别向各工地运送多少吨水泥，使总的吨千米数最小。
（2）为了进一步减少吨千米数，打算舍弃两个临时料场，改建两个新的，日储量各为20吨，问应建在何处，节省的吨千米数有多大？

• 问分析

记工地 $i$ 的坐标为$(a_i,b_i),i=1,2,...,6$ 水泥日用量为$d_i$
记料场 $j$ 的坐标为 $(x_j,y_j),j=1,2$ 日储量为20
记从料场 $j$ 往工地 $i$ 运送的水泥为 $x_{ij}$ 吨
则
$$\begin{gathered} \min z={\sum }_{i=1}^6{\sum }_{j=1}^2{x}_{ij}\sqrt{(a_i-x_i{)}^2+(b_i-y_i{)}^2} \\ s.t.\{\begin{array}{ll}{\sum }_{j=1}^2{x}_{ij}=d_i& i=1,2,...,6\\ {\sum }_{i=1}^6{x}_{ij}\le 20& j=1,2\\ x_{ij}\ge 0& i=1,2,..,6;j=1,2\end{array} \end{gathered}$$
（1）问$(x_j,y_j),j=1,2$是已知的，故共有12个变量，（2）问有16个变量
决策变量双下下标改为单下表：

$x_1$	$x_2$	…	$x_6$	$x_7$	$x_8$	…	$x_{12}$	$x_{13}$	$x_{14}$	$x_{15}$	$x_{16}$
$x_{11}$	$x_{21}$	…	$x_{61}$	$x_{12}$	$x_{22}$	…	$x_{62}$	$x_1$	$y_1$	$x_2$	$y_2$

function f = fun5(xx)  % 注意为了避免和下面的x同号，我们把决策变量的向量符号用xx表示（注意xx的长度为16）
    a=[1.25  8.75  0.5  5.75  3  7.25];  % 工地的横坐标
    b=[1.25  0.75  4.75	5  6.5  7.25];   % 工地的纵坐标
    x = [xx(13)  xx(15)];  % 新料场的横坐标
    y = [xx(14)  xx(16)];  % 新料场的纵坐标
    c = [];  % 初始化用来保存工地和料场距离的向量 (这个向量就是我们的系数向量）
    for  j =1:2
        for i = 1:6
            c = [c;  sqrt( (a(i)-x(j))^2 + (b(i)-y(j))^2)];  % 每循环一次就在c的末尾插入新的元素
        end
    end
    % 下面我们要求吨千米数，注意c是列向量，我们计算非线性规划时给定的初始值x0是行向量
    f = xx(1:12) * c;
end

clear;clc
format long g
% (1) 不等式约束
A =zeros(2,16);  % 注意这里要改成16
A(1,1:6) = 1;
A(2,7:12) = 1;
b = [20,20]';
% (2) 等式约束
Aeq = zeros(6,16);  % 注意这里要改成16
for i = 1:6
    Aeq(i,i) = 1;  Aeq(i,i+6) = 1;
end
beq = [3 5 4 7 6 11]';  % 每个工地的日需求量
%（3）上下界
lb = zeros(16,1);
% lb = [zeros(12,1); -inf*ones(4,1)];  两个新料场坐标的下界可以设为-inf

% 进行求解
x0 = [3 5 0 7 0 1 0 0 4 0 6 10 5 1 2 7];  % 用第一问的结果作为初始值
[x,fval] = fmincon(@fun5,x0,A,b,Aeq,beq,lb)  % 注意没有非线性约束，所以这里可以用[]替代，或者干脆不写

x =

1 至 4 列

      2.93079373691722          4.62098894690306          3.86630967073392           6.9323464081422

5 至 8 列

      1.54452626271135         0.025293774364618        0.0692062630827831         0.379011053096942

9 至 12 列

     0.133690329266078        0.0676535918578045          4.45547373728865          10.9747062256354

13 至 16 列

      5.74190837844137          4.99144031310414          7.24999993630655          7.24999999514472

fval =

      89.9232875546202