前言
通过模型算法，熟练对Matlab和python的应用。
学习视频链接：
https://www.bilibili.com/video/BV1EK41187QF?p=24&vd_source=67471d3a1b4f517b7a7964093e62f7e6
一、非线性规划法
- 非线性规划是一种求解目标函数或约束条件中有一个或几个非线性函数的最优化问题的方法运筹学的一个重要分支。20世纪50年代初，库哈（H.W.Kuhn) 和托克 (A.W. Tucker）提出了非线性规划的基本定理，为非线性规划奠定了理论基础。
- 特点：模型中至少一个变量是非线性
- 非线性规划模型的标准型：
m i n f ( x ) s.t. { A x ≤ b , A e q ⋅ x = b e q (线性) c ( x ) ≤ 0 , C e q ( x ) = 0 (非线性) l b ≤ x ≤ u b min\quad f(x)\\[1em]\\\text{s.t.}\begin{cases}Ax\leq b, Aeq\cdot x=beq&\text{(线性)}\\c(x)\leq0, Ceq(x)=0&\text{(非线性)}\\lb\leq x\leq ub\end{cases} minf(x)s.t.⎩ ⎨ ⎧Ax≤b,Aeq⋅x=beqc(x)≤0,Ceq(x)=0lb≤x≤ub(线性)(非线性)
二、例题：选址问题
- 临时料场：A(5,1)，B(2,7)；日储量各20吨
- (1）试制定每天的供应计划，即从两料场分别向各工地运送多少吨水泥，使总的吨千米数最小？
- (2）为了进一步减少吨千米数，打算舍弃两个临时料场，改建两个新的，日储量各为20吨，问应建在何处，节省的吨千米数为多大？（吨千米数：吨 * 千米）
1.确定决策变量
- 设第 i 个工地的坐标 ( a i , b i ) (a_i,b_i) (ai,bi) ，水泥日用量 d i d_i di ， i = 1 , 2 , . . . , 6 i=1,2,...,6 i=1,2,...,6 ；料场位置 ( x i ， y i ) (x_i，y_i) (xi，yi)，日储量 e j e_j ej， j = 1 , 2 j=1,2 j=1,2；从料场 j 向工地 i 的运送量为 x i j x_{ij} xij
2.确定约束条件
- 料场水泥运输总量不超过其日储量： ∑ i = 1 6 x i j ≤ e j , j = 1 , 2 \sum_{i=1}^{6}x_{ij}\leq e_{j},j=1,2 ∑i=16xij≤ej,j=1,2
- 两个料场向某工地运输量之和等于该工地水泥日用量： ∑ j = 1 2 x i j = d i , i = 1 , 2 , ⋯ , 6 \sum_{j=1}^{2}x_{ij}=d_{i},i=1,2,\cdots,6 ∑j=12xij=di,i=1,2,⋯,6
3.确定目标函数
- 求总吨千米数最小，即运送量乘运送距离求和最小 m i n f = ∑ j = 1 2 ∑ i = 1 6 x i j ( x j − a i ) 2 + ( y j − b i ) 2 min\quad f=\sum_{j=1}^{2}\sum_{i=1}^{6}x_{ij}\sqrt{\left(x_{j}-a_{i}\right)^{2}+\left(y_{j}-b_{i}\right)^{2}} minf=j=1∑2i=1∑6xij(xj−ai)2+(yj−bi)2
4.建立模型
m i n f = ∑ j = 1 2 ∑ i = 1 6 x i j ( x j − a i ) 2 + ( y j − b i ) 2 s . t . { ∑ i = 1 6 x i j ≤ e j , j = 1 , 2 ∑ j = 1 2 x i j = d i , i = 1 , 2 , … , 6 x i j ≥ 0 , i = 1 , 2 , … , 6 ; j = 1 , 2 \begin{aligned}&min\quad f=\sum_{j=1}^2\sum_{i=1}^6x_{ij}\sqrt{\left(x_j-a_i\right)^2+\left(y_j-b_i\right)^2}\\[0em]\\&{s}.t.\begin{cases}\sum_{i=1}^6x_{ij}\leq e_j,j=1,2\\[0em]\\\sum_{j=1}^2x_{ij}=d_i,i=1,2,\ldots,6\\[0em]\\x_{ij}\geq0,i=1,2,\ldots,6;j=1,2\end{cases}\end{aligned} minf=j=1∑2i=1∑6xij(xj−ai)2+(yj−bi)2s.t.⎩ ⎨ ⎧∑i=16xij≤ej,j=1,2∑j=12xij=di,i=1,2,…,6xij≥0,i=1,2,…,6;j=1,2
5.求解
- 对于第一问：因料场位置已知，故决策变量仅为 x i j x_{ij} xij ，为线性规划模型
- 对于第二问：新料场位置未知，所以 x j x_j xj 和 y j y_j yj 均为变量，且不是线性的，故为非线性规划模型
- 共有 8 个约束
m i n f = ∑ j = 1 2 ∑ i = 1 6 x i j ( x j − a i ) 2 + ( y j − b i ) 2 s t . { ∑ i = 1 6 x i j ≤ e j , j = 1 , 2 ( x 11 + x 21 + … + x 61 ≤ e 1 , x 12 + x 22 + … + x 62 ≤ e 2 ) ∑ j = 1 2 x i j = d i , i = 1 , 2 , … , 6 ( x 11 + x 12 = d 1 , … , x 61 + x 62 = d 6 ) x i j ≥ 0 , i = 1 , 2 , … , 6 ; j = 1 , 2 min \quad f=\sum_{j=1}^{2}\sum_{i=1}^{6}x_{ij}\sqrt{\left(x_{j}-a_{i}\right)^{2}+\left(y_{j}-b_{i}\right)^{2}}\\[1em]\\\\st.\begin{cases}\sum_{i=1}^{6}x_{ij}\leq e_{j},j=1,2 \quad\left(x_{11}+x_{21}+\ldots+x_{61}\leq e_{1},x_{12}+x_{22}+\ldots+x_{62}\leq e_{2}\right)\\\sum_{j=1}^{2}x_{ij}=d_{i},i=1,2,\ldots,6\quad\left(x_{11}+x_{12}=d_{1},\ldots,x_{61}+x_{62}=d_{6}\right)\\x_{ij}\geq0,i=1,2,\ldots,6;j=1,2\end{cases} minf=j=1∑2i=1∑6xij(xj−ai)2+(yj−bi)2st.⎩ ⎨ ⎧∑i=16xij≤ej,j=1,2(x11+x21+…+x61≤e1,x12+x22+…+x62≤e2)∑j=12xij=di,i=1,2,…,6(x11+x12=d1,…,x61+x62=d6)xij≥0,i=1,2,…,6;j=1,2
三、代码实现----Matlab
1.Matlab 的 fmincon 函数
（1）基本用法
fmincon 的基本调用格式如下：
[x, fval] = fmincon(fun, x0, A, b, Aeq, beq, lb, ub, nonlcon, options)
输入参数
- fun: 目标函数的函数句柄。把目标函数定义为一个单独的函数文件(min)
- x0: 初始点，决策变量的初始猜测值。
- A, b: 线性不等式约束，形式为
A*x <= b（标准型为 ≤ ）。 - Aeq, beq: 线性等式约束，形式为
Aeq*x = beq。 - lb, ub: 决策变量的下界和上界。
- nonlcon: 非线性约束的函数句柄。例如，
@mycon，其中mycon是一个定义非线性约束的 MATLAB 函数，应返回两个输出：c和ceq，分别表示非线性不等式和等式约束。 - options: 优化选项（默认内点法），通过
optimoptions函数创建。
输出参数
- x: 优化问题的解，决策变量的最优值。
- fval: 在最优解处的目标函数值。
fmincon 函数
[x, fval] = fmincon(fun, x0, A, b, Aeq, beq, lb, ub, nonlcon, options)
- 其中非线性规划中对于初始值x0的选取非常重要，因为非线性规划的算法求解出来的是一个局部最优解。如果要求全局最优解，有两个思路：
- 给定不同的初始值，在里面找到一个最优解；
- 先用蒙特卡罗模拟，得到一个蒙特卡罗解，然后将这个解作为初始值来求最优解
- option选项可以给定求解的算法，一共有五种，interior-point（内点法）、trust-region-ref lective（信赖域反射法）、sqp（序列二次规划法）、sqp-legacy（约束非线性优化算法）、acti ve-set（有效集法）。不同的算法有其各自的优缺点和适用情况 我们可以改变求解的算法来对比求解的结果。
- fun表示目标函数，我们要编写一个独立的 “m文件” 储存目标函数
function f = fun(x)
f = ......
end
- nonlcon 表示非线性部分的约束，也要编写一个独立的 “m文件” 储存非线性约束条件
function [c,ceq] = nonlfun(x)
c=［非线性不等式约束1;
.......
非线性不等式约束2]
ceq=［非线性等式约束1;
.......
非线性等式约束2]
end
- 若不存在某种约束，可以用 “[ ]” 替代，若后面全为 “[ ]” 且option使用默认，后面的 “[ ]” 可省略
（2）简单示例
m
i
n
y
=
x
1
2
+
x
2
2
−
x
1
x
2
−
2
x
1
−
5
x
2
,
s
.
t
.
{
−
(
x
1
−
1
)
2
+
x
2
≥
0
(
1
)
,
2
x
1
−
3
x
2
+
6
≥
0
(
2
)
\begin{gathered}min\quad {y}=x_1^2+x_2^2-x_1x_2-2x_1-5x_2,\\\mathrm{s.t.}\begin{cases}-\left(x_1-1\right)^2+x_2\geq0 \quad(1),\\2x_1-3x_2+6\geq0 \quad(2)&\end{cases}\end{gathered}
miny=x12+x22−x1x2−2x1−5x2,s.t.{−(x1−1)2+x2≥0(1),2x1−3x2+6≥0(2)
注：
标准型中是 ≤ ，所以不等式两边要同时乘 ‘-’
(1)式是非线性不等式；(2)式是线性不等式；
Matlab代码
clear;clc
format long g
% 可以将Matlab的计算结果显示为一般的长数字格式(默认会保留四位小数,或使用科学计数法)
x0 = [0,0];
A = [-2,3];
b = 6;
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlcon1);
%%
% 使用SQP算法(序列二次规划)
clear;clc;
x0 = [0,0];
A = [-2,3];
b = 6;
option = optimoptions('fmincon','Algorithm','sqp');
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlcon1);
%%
% 使用蒙特卡洛的方法来找初始值
clear;clc
format long g
n = 100000;
x1 = unifrnd(-100,100,n,1); % 生成在[-100,100]范围内均匀分布的随机数组组成的n行1列的向量
x2 = unifrnd(-100,100,n,1);
y0 = +inf;
for i = 1:n
if (-(x1(i) - 1)^2 + x2(i)) >= 0 && (2*x1(i) - 3*x2(i) + 6) >=0
fval = x1(i)^2 + x2(i)^2 - x1(i)*x2(i) - 2*x1(i) - 5*x2(i);
if fval < y0
y0 = fval;
x = [x1(i), x2(i)];
end
end
end
disp('蒙特卡洛的方法找到初始值为')
disp(x)
x0 = x;
A = [-2,3];
b = 6;
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlcon1);
function f = fun1(x)
f = x(1)^2 + x(2)^2 - x(1) * x(2) - 2 * x(1) - 5 * x(2);
end
function [c,ceq] = nonlcon1(x)
% 注意:这里的c实际上就是非线性不等式约束,ceq实际上就是非线性等式约束
% 输入的x是决策变量
c = [(x(1) - 1)^2 - x(2)];
ceq = [];
end
运行结果：
默认算法（内点法）得到的结果：
使用SQP算法得到的结果：
加入蒙特卡洛法得到的结果：
2.Matlab 代码
第一问：线性规划 使用 linprog 函数
由于决策变量为 x i j x_{ij} xij ，所以变量有 12 个，将双角标改为单角标变量，组成一维列向量，即 x 11 → x 1 , x 21 → x 2 , … , x 62 → x 12 x_{11}\rightarrow x_{1} ,\quad x_{21}\rightarrow x_{2} ,\quad\ldots\quad,\quad x_{62}\rightarrow x_{12} x11→x1,x21→x2,…,x62→x12
% 第一问: 线性规划
clear;clc
% [x,fval] = linprog(f,A,b,Aeq,beq,lb,ub);
% f 目标函数的系数向量（必须是求最小值形式下的）
% A，b 不等式约束条件的变量系数矩阵和常数项矩阵（必须是 ≤ 形式）
% Aeq，beq 等式约束条件的系数矩阵和常数项矩阵
% lb，ub 决策变量的最小取值和最大取值
% x 是返回的最优解的变量取值， fval 返回目标函数的最优值
A = [[ones(1,6), zeros(1,6)];[zeros(1,6), ones(1,6)]];
b = [20;20];
Aeq = [eye(6),eye(6)];
beq = [3;5;4;7;6;11];
lb = zeros(12,1);
piont = repmat([1.25 1.25;
8.75 0.75;
0.5 4.75;
5.75 5;
3 6.5;
7.25 7.25],2,1);
position = [repmat([5 1],6,1);repmat([2 7],6,1)];
f = (sum((position - piont) .* (position - piont),2)) .^ 0.5;
[x,fval] = linprog(f,A,b,Aeq,beq,lb);
运行结果：
第二问：非线性规划
由于新料场位置未知，所以 x j x_j xj 和 y j y_j yj 均为变量，有两个新料场，并且每个料场坐有 x，y两个变量，所以相比较第一问，变量增加了四个，我将增加的四个变量加到了原本的 12 个变量之前，这样一维列向量的元素就变成 16 个。
clear;clc
%% 非线性规划的函数
% [x,fval] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
% x0 表示给定的初始值(用行向量或者列向量表示),必须得写
% A b 表示线性不等式约束
% Aeq beq 表示线性等式约束
% lb ub 表示上下界约束
% fun 表示目标函数
% nonlcon 表示非线性约束的函数
% option 表示求解非线性规划使用的方法
% 初始值是第一问求得的数值，前四个变量也是第一问题目中的数值
x0 = [5;1;2;7;3;5;0;7;0;1;0;0;4;0;6;10];
A = [zeros(2,4),[[ones(1,6), zeros(1,6)];[zeros(1,6), ones(1,6)]]];
b = [20;20];
Aeq = [zeros(6,4),[eye(6),eye(6)]];
beq = [3;5;4;7;6;11];
lb = [-inf; -inf; -inf; -inf; zeros(12,1)];
[x,fval] = fmincon(@fun2,x0,A,b,Aeq,beq,lb);
function ff = fun2(x) % x1:x(1) y1:x(2) x2:x(3) y2:x(4)
piont = repmat([1.25 1.25;
8.75 0.75;
0.5 4.75;
5.75 5;
3 6.5;
7.25 7.25],2,1);
position = [repmat([x(1) x(2)],6,1);repmat([x(3) x(4)],6,1)];
ff = sum((sum((position - piont) .* (position - piont),2)) .^ 0.5 .* x(5:end,1));
end
运行结果：
第二问：非线性规划(加入蒙特卡洛法找初始值）
% 第二问使用蒙特卡洛求近似解
% 1.因为第2问模型中有16个变量，所以要给16个变量分别生成随机值，作为当前解；
% 2.判断这些当前解是否满足模型的约束条件
% 3.若满足，代入目标函数，求当前目标函数值
% 4.判断当前目标函数值是否比已求的较优函数值更好，若是，则替换掉较优函数值和对应的较优解
% 5.不断重复前4步，构成统计意义，求得较优解。
%%
% 后12个数是6个工地从两个料场接收的量，6个工地分别需要3,5,4,7,6,11吨水泥
% 所以后12个变量分别需要取0到3，0到5，……，0到11的随机数
% 为加速求近似，取后12个变量（工地从料场接受的量）为随机整数
% randi(n)为随机取1到n之间的一个整数，则randi(n+1)-1为取0到n间随机整数
% 前4个变量是两个料场的横纵坐标，取一定范围内的随机数（根据题目，工地坐标都在0到9，所以此处取0到9）
p0 = 10000;
n = 10^6;
x_m0 = 0;
c = [zeros(6,4),[eye(6),eye(6)]];
c_1 = [3;5;4;7;6;11];
for i = 1:n
x = [randi(10)-1; randi(10)-1; randi(10)-1; randi(10)-1; ...
randi(4)-1; randi(6)-1; randi(5)-1; randi(8)-1; randi(7)-1; randi(12)-1;...
randi(4)-1; randi(6)-1; randi(5)-1; randi(8)-1; randi(7)-1; randi(12)-1];
cd1 = sum(x(5:10,:));
cd2 = sum(x(11:16,:));
cd3 = c * x - c_1; % 按列求和得到列向量 6*1
if cd1 <= 20 && cd2 <= 20 && all(abs(cd3)<=1)
ff = fun2(x);
if ff < p0
p0 = ff;
x_m0 = x;
end
end
end
res1 = ff;
res2 = x_m0;
A = [zeros(2,4),[[ones(1,6), zeros(1,6)];[zeros(1,6), ones(1,6)]]];
b = [20;20];
Aeq = [zeros(6,4),[eye(6),eye(6)]];
beq = [3;5;4;7;6;11];
lb = [-inf; -inf; -inf; -inf; zeros(12,1)];
[x,fval] = fmincon(@fun2,x_m0,A,b,Aeq,beq,lb);
运行结果：
可以发现使用蒙特卡洛法找到初始值后，得到的吨千米数比直接调用fmincon 小，即加入蒙特卡洛法效果更好。
四、代码实现----python
1.python 的 minimize 函数
在 MATLAB 中，fmincon 是一个非常常用的函数，用于求解带约束的非线性优化问题。在 Python 中，类似功能可以通过 SciPy 库中的 scipy.optimize.minimize 函数实现。SciPy 提供了多种优化算法，并且支持处理等式和不等式约束、边界条件等。
scipy.optimize.minimize 是 SciPy 库中的一个强大工具，用于求解无约束或带约束的最优化问题。它提供了多种最优化算法，可以处理各种类型的优化问题，包括非线性优化、带边界条件的优化和带约束条件的优化。下面是关于 minimize 函数的详细介绍和使用示例。
（1）基本用法
scipy.optimize.minimize 的基本调用格式如下：
scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None,
bounds=None, constraints=(), tol=None, callback=None, options=None)
参数
- fun: 需要最小化的目标函数。
- x0: 初始猜测值，是一个数组。
- args: 传递给目标函数和约束函数的额外参数。
- method: 指定使用的优化算法，如
'BFGS','CG','L-BFGS-B','TNC','SLSQP'等（默认'SLSQP'算法）。 - jac: 目标函数的梯度，可以是布尔值、数组或函数。
- hess: 目标函数的Hessian矩阵，用于二次优化算法。
- hessp: Hessian矩阵乘以向量，用于信赖域算法。
- bounds: 变量的边界条件，用于约束优化。
- constraints: 约束条件，可以是字典或字典列表。
- tol: 终止条件的公差。
- callback: 每次迭代后调用的函数。
- options: 其他选项，如最大迭代次数、显示信息等。
常用参数
- fun: 需要最小化的目标函数。
- x0: 初始猜测值，是一个数组。
- method: 指定使用的优化算法，如
'BFGS','CG','L-BFGS-B','TNC','SLSQP'等（默认'SLSQP'算法）。 - bounds: 变量的边界条件，用于约束优化。
- constraints: 约束条件，可以是字典或字典列表。
常用优化方法
- BFGS：Broyden–Fletcher–Goldfarb–Shanno 算法，用于无约束优化。
- L-BFGS-B：Limited-memory BFGS 算法，可以处理边界条件。
- TNC：Truncated Newton Conjugate-Gradient 算法，也能处理边界条件。
- SLSQP：Sequential Least SQuares Programming 算法，可以处理等式和不等式约束以及边界条件。
scipy.optimize.minimize函数返回一个OptimizeResult对象，该对象包含有关优化结果的各种信息。以下是OptimizeResult对象中包含的主要属性：
函数返回值
- x: 优化问题的解，即目标函数最小化时的变量值。
- success: 布尔值，表示优化是否成功。
- status: 整数，表示优化结束的状态码。
- message: 字符串，表示优化结束的详细信息。
- fun: 目标函数在解处的值。
- jac: 目标函数在解处的雅可比矩阵（梯度）。
- hess_inv: 目标函数的逆Hessian矩阵（如果可用）。
- nfev: 目标函数的评估次数。
- njev: 雅可比矩阵（梯度）的评估次数。
- nhev: Hessian矩阵的评估次数。
- nit: 迭代次数。
- maxcv: 最大约束违反值。
参考文档
更多关于 scipy.optimize.minimize 的详细信息和用法，可以参考 SciPy 官方文档。
（2）简单示例
m i n y = x 1 2 + x 2 2 − x 1 x 2 − 2 x 1 − 5 x 2 , s . t . { − ( x 1 − 1 ) 2 + x 2 ≥ 0 ( 1 ) , 2 x 1 − 3 x 2 + 6 ≥ 0 ( 2 ) \begin{gathered}min\quad {y}=x_1^2+x_2^2-x_1x_2-2x_1-5x_2,\\\mathrm{s.t.}\begin{cases}-\left(x_1-1\right)^2+x_2\geq0 \quad(1),\\2x_1-3x_2+6\geq0 \quad(2)&\end{cases}\end{gathered} miny=x12+x22−x1x2−2x1−5x2,s.t.{−(x1−1)2+x2≥0(1),2x1−3x2+6≥0(2)
值得注意的是，与 Matlab 的fmincon函数不同，scipy.optimize.minimize 中不等式约束函数目标是大于等于0的，所以不需要乘负号。
python代码
import numpy as np
from scipy.optimize import minimize
# 定义目标函数
def fun(x):
return x[0]**2 + x[1]**2 - x[0]*x[1] - 2*x[0] - 5*x[1]
# 初始猜测数组
x0 = np.array([0,0])
# eq表示 本约束函数目标等于0 ； ineq 表示 约束函数目标大于等于0
# 不等式约束1
def ineq_constraint_1(x):
return -(x[0]-1)**2 + x[1]
# 不等式约束2
def ineq_constraint_2(x):
return 2*x[0] - 3*x[1] + 6
# 定义约束
constraints = [
{'type': 'ineq', 'fun': ineq_constraint_1}, # 不等式约束1
{'type': 'ineq', 'fun': ineq_constraint_2} # 不等式约束2
]
result = minimize(fun=fun, x0=x0, constraints=constraints)
print(result.x)
print(result.fun)
运行结果：
# 使用蒙特卡洛寻找初始值
import numpy as np
from scipy.optimize import minimize
n = 100000
x1 = np.random.uniform(-100, 100, (n, 1))
x2 = np.random.uniform(-100, 100, (n, 1))
y = +np.inf
for i in range(n):
if -(x1[i]-1)**2 + x2[i] >= 0 and 2*x1[i] - 3*x2[i] + 6 >= 0:
p = x1[i]**2 + x2[i]**2 - x1[i]*x2[i] - 2*x1[i] - 5*x2[i]
if p < y:
y = p
x3 = [float(x1[i]),float(x2[i])]
print(x3)
def fun(x):
return x[0]**2 + x[1]**2 - x[0]*x[1] - 2*x[0] - 5*x[1]
# eq表示 本约束函数目标等于0 ； ineq 表示 约束函数目标大于等于0
# 不等式约束1
def ineq_constraint_1(x):
return -(x[0]-1)**2 + x[1]
# 不等式约束2
def ineq_constraint_2(x):
return 2*x[0] - 3*x[1] + 6
# 定义约束
constraints = [
{'type': 'ineq', 'fun': ineq_constraint_1}, # 不等式约束1
{'type': 'ineq', 'fun': ineq_constraint_2} # 不等式约束2
]
result = minimize(fun=fun, x0=x3, constraints=constraints)
print(result.x)
print(result.fun)
运行结果：
2.python 代码
第一问：线性规划 使用 linprog 函数
import numpy as np
from scipy.optimize import linprog
# 使用 np.hstack 将两个1x6的数组水平堆叠,形成一个1x12的数组。
# 使用 np.vstack 将两个1x12的数组垂直堆叠,形成一个2x12的数组。
A = np.vstack([np.hstack([np.ones((1,6)), np.zeros((1,6))]),np.hstack([np.zeros((1,6)), np.ones((1,6))])])
b = np.array([[20], [20]])
A_eq = np.hstack([np.eye(6),np.eye(6)])
b_eq = np.array([[3], [5], [4], [7], [6], [11]])
bounds = [(0,None)] * 12
piont = np.tile(np.array([[1.25,1.25],
[8.75,0.75],
[0.5,4.75],
[5.75,5],
[3,6.5],
[7.25,7.25]]),(2,1))
position = np.vstack([np.tile(np.array([5,1]),(6,1)),np.tile(np.array([2,7]),(6,1))])
f = (np.sum((position - piont) * (position - piont),1)) ** 0.5
res = linprog(f,A,b,A_eq,b_eq,bounds)
print(res.x)
print(res.fun)
运行结果：
第二问：非线性规划
import numpy as np
from scipy.optimize import minimize
x0 = np.array([5, 1, 2, 7, 3, 5, 0, 7, 0, 1, 0, 0, 4, 0, 6, 10])
# 不等式约束
def ineq_constraint(x):
A = np.hstack([np.zeros((2, 4)),
np.vstack([np.hstack([np.ones((1, 6)),
np.zeros((1, 6))]),
np.hstack([np.zeros((1, 6)),
np.ones((1, 6))])])])
b = np.array([20, 20]) # 修改为一维数组
return - (A @ x).flatten() + b # 展平为一维数组
# 等式约束
def eq_constraint(x):
A_eq = np.hstack([np.zeros((6, 4)), np.hstack([np.eye(6), np.eye(6)])])
b_eq = np.array([3, 5, 4, 7, 6, 11]) # 修改为一维数组
return (A_eq @ x).flatten() - b_eq # 展平为一维数组
# 定义约束
constraints = [
{'type': 'eq', 'fun': eq_constraint}, # 等式约束
{'type': 'ineq', 'fun': ineq_constraint} # 不等式约束
]
bounds = [(None, None)] * 4 + [(0, None)] * 12
def fun(x):
piont = np.tile(np.array([[1.25, 1.25],
[8.75, 0.75],
[0.5, 4.75],
[5.75, 5],
[3, 6.5],
[7.25, 7.25]]), (2, 1))
position = np.vstack([np.tile(np.array([x[0], x[1]]), (6, 1)),
np.tile(np.array([x[2], x[3]]), (6, 1))])
ff = np.sum((np.sum((position - piont) * (position - piont), 1)) ** 0.5 * x[4:])
return ff
result = minimize(fun=fun, x0=x0, bounds=bounds, constraints=constraints)
print(result.x)
print(result.fun)
运行结果：
第二问：非线性规划(加入蒙特卡洛法找初始值）
import numpy.random as rd
import numpy as np
from scipy.optimize import minimize
# 不等式约束
def ineq_constraint(x):
A = np.hstack([np.zeros((2, 4)),
np.vstack([np.hstack([np.ones((1, 6)),
np.zeros((1, 6))]),
np.hstack([np.zeros((1, 6)),
np.ones((1, 6))])])])
b = np.array([20, 20]) # 修改为一维数组
return - (A @ x).flatten() + b # 展平为一维数组
# 等式约束
def eq_constraint(x):
A_eq = np.hstack([np.zeros((6, 4)), np.hstack([np.eye(6), np.eye(6)])])
b_eq = np.array([3, 5, 4, 7, 6, 11]) # 修改为一维数组
return (A_eq @ x).flatten() - b_eq # 展平为一维数组
# 定义约束
constraints = [
{'type': 'eq', 'fun': eq_constraint}, # 等式约束
{'type': 'ineq', 'fun': ineq_constraint} # 不等式约束
]
bounds = [(None, None)] * 4 + [(0, None)] * 12
def fun(x):
piont = np.tile(np.array([[1.25, 1.25],
[8.75, 0.75],
[0.5, 4.75],
[5.75, 5],
[3, 6.5],
[7.25, 7.25]]), (2, 1))
position = np.vstack([np.tile(np.array([x[0], x[1]]), (6, 1)),
np.tile(np.array([x[2], x[3]]), (6, 1))])
ff = np.sum((np.sum((position - piont) * (position - piont), 1)) ** 0.5 * x[4:])
return ff
p0 = 10000
n = 10 ** 6
x_m0 = np.array(np.zeros((16,)))
c = np.hstack([np.zeros((6,4)),np.hstack([np.eye(6),np.eye(6)])]) # 6*16
c_1 = np.array([3, 5, 4, 7, 6, 11])
for i in range(0,n):
x = np.array([rd.randint(0,9), rd.randint(0,9), rd.randint(0,9), rd.randint(0,9),
rd.randint(0,3), rd.randint(0,5), rd.randint(0,4), rd.randint(0,7), rd.randint(0,6), rd.randint(0,11),
rd.randint(0,3), rd.randint(0,5), rd.randint(0,4), rd.randint(0,7), rd.randint(0,6), rd.randint(0,11)])
cd1 = np.sum(x[4:10])
cd2 = np.sum(x[10:16])
cd3 = c @ x - c_1 # 按列求和得到列向量 6*1
if cd1 <= 20 and cd2 <= 20 and np.all(np.abs(cd3) <= 1):
ff = fun(x)
if ff < p0:
p0 = ff
x_m0 = x
res1 = x_m0
result = minimize(fun=fun, x0=res1, bounds=bounds, constraints=constraints)
print(result.x)
print(result.fun)
运行结果：
比较两种方法同样得出结论：加入蒙特卡洛法效果更好。
遇到的问题
错误代码：
import numpy as np
from scipy.optimize import minimize
x0 = np.array([5, 1, 2, 7, 3, 5, 0, 7, 0, 1, 0, 0, 4, 0, 6, 10])
# eq表示 本约束函数目标等于0 ； ineq 表示 约束函数目标大于等于0
# 不等式约束
def ineq_constraint(x):
A = np.hstack([np.zeros((2,4)),
np.vstack([np.hstack([np.ones((1,6)),
np.zeros((1,6))]),
np.hstack([np.zeros((1,6)),
np.ones((1,6))])])])
b = np.array([[20], [20]])
return - (A @ x).reshape(-1,1) - b
# 等式约束
def eq_constraint(x):
A_eq = np.hstack([np.zeros((6,4)),np.hstack([np.eye(6),np.eye(6)])])
b_eq = np.array([[3], [5], [4], [7], [6], [11]])
return (A_eq @ x).reshape(-1,1) - b_eq
# 定义约束
constraints = [
{'type': 'eq', 'fun': eq_constraint}, # 等式约束
{'type': 'ineq', 'fun': ineq_constraint} # 不等式约束
]
bounds = [(None,None)] * 4 + [(0,None)] * 12
def fun(x):
piont = np.tile(np.array([[1.25, 1.25],
[8.75, 0.75],
[0.5, 4.75],
[5.75, 5],
[3, 6.5],
[7.25, 7.25]]),(2,1))
position = np.vstack([np.tile(np.array([x[0],x[1]]),(6,1)),np.tile(np.array([x[2],x[3]]),(6,1))])
ff = np.sum((np.sum((position - piont) * (position - piont),1)) ** 0.5 * x[4:])
return ff
result = minimize(fun=fun, x0=x0, bounds=bounds, constraints=constraints)
错误之处是在定义 不等式约束 和 等式约束时使用了 reshape(-1,1)。在未对数组进行reshape(-1,1)时，报错：
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[1], line 46
43 ff = np.sum((np.sum((position - piont) * (position - piont),1)) ** 0.5 * x[4:])
44 return ff
---> 46 result = minimize(fun=fun, x0=x0, bounds=bounds, constraints=constraints)
File c:\Users\zhy\AppData\Local\Programs\Python\Python39\lib\site-packages\scipy\optimize\_minimize.py:722, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
719 res = _minimize_cobyla(fun, x0, args, constraints, callback=callback,
720 bounds=bounds, **options)
721 elif meth == 'slsqp':
--> 722 res = _minimize_slsqp(fun, x0, args, jac, bounds,
723 constraints, callback=callback, **options)
724 elif meth == 'trust-constr':
725 res = _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
726 bounds, constraints,
727 callback=callback, **options)
File c:\Users\zhy\AppData\Local\Programs\Python\Python39\lib\site-packages\scipy\optimize\_slsqp_py.py:426, in _minimize_slsqp(func, x0, args, jac, bounds, constraints, maxiter, ftol, iprint, disp, eps, callback, finite_diff_rel_step, **unknown_options)
424 fx = wrapped_fun(x)
425 g = append(wrapped_grad(x), 0.0)
--> 426 c = _eval_constraint(x, cons)
427 a = _eval_con_normals(x, cons, la, n, m, meq, mieq)
429 while 1:
430 # Call SLSQP
...
487 # Now combine c_eq and c_ieq into a single matrix
--> 488 c = concatenate((c_eq, c_ieq))
489 return c
ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 6 and the array at index 1 has size 2
报错原因是：在调用minimize函数时，处理等式约束和不等式约束得到的结果数组，在尝试拼接两个数组时，这两个数组在拼接轴以外的其他轴上的尺寸不一致。
# Now combine c_eq and c_ieq into a single matrix
c = concatenate((c_eq, c_ieq))
return c
但加入 reshape(-1,1)后，报了另一个错误：
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
Cell In[2], line 46
43 ff = np.sum((np.sum((position - piont) * (position - piont),1)) ** 0.5 * x[4:])
44 return ff
---> 46 result = minimize(fun=fun, x0=x0, bounds=bounds, constraints=constraints)
File c:\Users\zhy\AppData\Local\Programs\Python\Python39\lib\site-packages\scipy\optimize\_minimize.py:722, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
719 res = _minimize_cobyla(fun, x0, args, constraints, callback=callback,
720 bounds=bounds, **options)
721 elif meth == 'slsqp':
--> 722 res = _minimize_slsqp(fun, x0, args, jac, bounds,
723 constraints, callback=callback, **options)
724 elif meth == 'trust-constr':
725 res = _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
726 bounds, constraints,
727 callback=callback, **options)
File c:\Users\zhy\AppData\Local\Programs\Python\Python39\lib\site-packages\scipy\optimize\_slsqp_py.py:427, in _minimize_slsqp(func, x0, args, jac, bounds, constraints, maxiter, ftol, iprint, disp, eps, callback, finite_diff_rel_step, **unknown_options)
425 g = append(wrapped_grad(x), 0.0)
426 c = _eval_constraint(x, cons)
--> 427 a = _eval_con_normals(x, cons, la, n, m, meq, mieq)
429 while 1:
430 # Call SLSQP
431 slsqp(m, meq, x, xl, xu, fx, c, g, a, acc, majiter, mode, w, jw,
...
--> 472 raise RuntimeError("`fun` return value has "
473 "more than 1 dimension.")
474 return f
RuntimeError: `fun` return value has more than 1 dimension.
报错原因是：目标函数fun的返回值具有多于1个维度。scipy.optimize.minimize函数要求目标函数返回一个标量值，而不是一个数组或矩阵。
定义 不等式约束 和 等式约束时加入 reshape(-1,1)后，数组就变成了二维数组，所以代码不能继续运行下去。
解决办法：使用 .flatten() 或 .ravel() 方法将数组展平成一维。
总结
本文介绍了非线性规划模型，并比较了加入蒙特卡洛法与不加入的区别，分别使用Matlab和python进行代码编写。
3388
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
