【学习笔记】Matlab和python双语言的学习(多目标规划）

文章目录

前言
一、多目标规划
- 1.特点
- 2.一般形式
- 3.多目标规划的解
- 4.线性加权法
二、典型示例-----化工厂生产问题
三、代码实现----Matlab
四、代码实现----python
总结

前言

通过模型算法，熟练对Matlab和python的应用。
学习视频链接：
https://www.bilibili.com/video/BV1EK41187QF?p=30&vd_source=67471d3a1b4f517b7a7964093e62f7e6

一、多目标规划

1.特点

衡量一个方案的角度往往不止一个

经济、管理、军事、科学和工程设计等领域衡量一个方案的好坏难以指标来判断，而需要用多个目标来比较，而这些目标有时不甚协调，甚至是矛盾的。
设计一个导弹，既要射程远，命中率高，还要耗燃料少
选择新厂址，除了要考虑运费、造价、燃料供应费等经济指标外，还要考虑对环境的污染社会因素
多目标规划是数学规划的一个分支。研究多于一个的目标函数在给定区域上的最优化。又称多目标最优化。通常记为 MQP(multi-objective programming)。
多目标规划的概念是 1961年由美国数学家查尔斯和库柏首先提出的。多目标最优化思想，最早是在1896年由法国经济学家V．帕雷托提出来的。他从政治经济学的角度考虑把本质上是不可比较的许多目标化成单个目标的最优化问题，从而涉及了多目标规划问题和多目标的概念

2.一般形式

多目标规划是多目标决策的重要内容之一，在进行多目标决策时，当希望每个目标都尽可能的大（或尽可能的小）时，就形成了一个多目标规划问题，其一般形式为：
$\begin{array}{rl}{\mathrm{min~f(x)=[f_1(x),~f_2(x),\cdots,f_m(x)]^T,}}\\{\mathrm{s.t.~\begin{cases}g_i(x)\leq0,&\mathrm{i=1,2,\cdots,p,}\\h_j(x)=0,&\mathrm{j=1,2,\cdots,q,}\end{cases}}}\end{array}$

其中 x 为决策向量， $f_1(x)$ ， $f_2(x)$ ，… ， $f_m(x)$ 为目标函数， $s . t .$ 式为约束条件
$记\quad\Omega=\left\{x|g_{i}(x)\leq0, i=1,2,\cdots,p; h_{j}(x)=0, j=1,2,\cdots,q\right. \}$
称 $\Omega$ 为多目标规划的可行域（决策空间）， $f(\Omega)=\left\{f(x)|x\in\Omega\right. \}$ 为多目标规划问题的像集（目标空间），多目标规划问题以下简称问题（MP)。

3.多目标规划的解

最优解：

帕累托最优解（Pareto Optimal Solution）描述了在给定多个目标函数的情况下，一个解的状态，使得没有其他解能够在不使至少一个目标变得更差的情况下使得一个或多个目标变得更好。
在多目标优化问题中，我们通常有多个目标函数 $f_1(x),f_2(x),\ldots,f_m(x)$ ，需要同时优化这些目标函数。假设 ( x ) 是决策变量的向量，定义如下：
一个解 $x^*$ 被称为帕累托最优解，如果不存在其他解 $x$ 使得以下两点同时成立：
1. 对于所有目标函数 $i$ ，有 $f_i(x)\leq f_i(x^*)$ ；
2. 至少对于一个目标函数 $j$ ，有 $f_j(x)<f_j(x^*)$
换句话说，一个解 $x^*$ 是帕累托最优的，如果没有其他解在所有目标函数上都优于 $x^*$ 且在至少一个目标函数上严格优于 $x^*$ 。

有效解：

帕累托有效解（Pareto Efficient Solution）在多目标优化中是一个非常重要的概念。它描述了在优化过程中那些不可被其他解支配的解，也就是说，在所有目标函数上没有其他解优于它们。
在多目标优化中，有多个目标函数 $f_1(x),f_2(x),\ldots,f_m(x)$ 需要同时优化。一个解 $x^*$ 被称为帕累托有效解，如果不存在其他解 $x$ 使得：
1. $f_i(x)\leq f_i(x^*)$ 对所有 $\ldots, m$
2. 且至少存在一个 $j$ 使得 $f_j(x) < f_j(x^*)$
换句话说，一个解 $x^*$ 是帕累托有效的，当且仅当没有其他解在所有目标函数上都优于 $x^*$ ，且在至少一个目标函数上严格优于 $x^*$ 。

满意解：

主要是从决策过程角度,根据决策者的偏好与要求而提出的，设可行域为 $\Omega$ ，要求 $m$ 个目标函数 $f_i(i=1,2,\cdots,m)$ 越小越好，有时决策者的期望较低,给出了 $m$ 个阈值 $\alpha_i$ ，当 $\bar{x}\in\Omega$ 满足 $f_i(\bar{x})\leq\alpha_i(i=1,2,\cdots,m)$ 时，就认为 $\bar{x}$ 是可以接受的、是满意的。这样的 $\bar{x}$ 就称为一个满意解。

多目标规划的求解

值得注意的是，在多目标规划中，除去目标函数一般是彼此冲突外，还有另一个特点：目标函数的不可公度性。所以通常在求解前，先对目标函数进行预处理。预处理的内容包括：无量纲化处理，归一化处理等

求解问题（MP）的有效解的常用方法

在这里插入图片描述

4.线性加权法

若一个规划问题中有多个目标，我们可以对多目标函数进行加权组合，使问题变为单目标规划，然后再利用之前学习的解法进行求解
线性加权法具体步骤如下：

第一步：确定每个目标函数的权系数。
$0\leq w_{j}\leq1, j=1,2,\cdots,m;\quad\sum_{j=1}^{m}w_{j}=1.$
第二步：写出评价函数 $\sum_{j=1}^{m}w_{j}f_{j}=1$
第三步：求评价函数最优值

$\sum_{i=1}^{m}w_{i}f_{i}(x),$
$\mathrm{s.t.}\quad x\in{\Omega}.$

二、典型示例-----化工厂生产问题

某化工厂今年拟生产两种新产品 A 和 B ，其生产费用分别为 2 万元／吨和 5 万元／吨。这两种产品均将造成环境污染，每生产一吨 A 产品会产生 0.4 吨的污染，每生产一吨 B 产品会产生 0. 3 吨的污染。由于条件限制，工厂生产产品 A 和 B 的最大生产能力各为每月 5 吨和 6 吨，而市场需要这两种产品的总量每月不少于 7 吨。该工厂决策认为，这两个目标中环境污染应该优先考虑，且根据经验生产费用的参考值为 30 万元，污染量参考值为 2 吨。试问工厂如何安排生产计划，在满足市场需要的前提下，使设备的花费和产生的污染均达到最小。

解：设工厂每月产品A生产 $x_1$ 吨，B生产 $x_2$ 吨，那么产生的污染分别为 $0. 4x_1$ 吨和 $0.3x_2$ 吨

建立多目标规划模型：
$\begin{cases}\min f_1=2x_1+5x_2\\\min f_2=0.4x_1+0.3x_2\\st.\begin{cases}x_1+{x}_2\geq7\\0\leq x_1\leq5\\0\leq x_2\leq6\end{cases}\end{cases}$

下面我们将其转换为一个单目标规划问题，即对上面的两个目标函数进行加权，由于该工厂决策认为环境污染应优先考虑，因此我们可以选取 $f_1$ 和 $f_2$ 的权重分别为 0.4 和 0.6 。注意到两个目标函数的单位不同，一个为"万元"，一个为"吨"，因此我们需要首先对目标函数进行标准化来消除量纲的影响，然后再进行加权，由于题目中已经给了产品费用和污染量的参考值，因此我们将这两个目标函数分别除以其参考值来消除量纲。

加权组合后的目标函数：
$f=0.4\times\frac{f_1}{30}+0.6\times\frac{f_2}{2}=\frac{0.4}{30}\times(2x_1+5x_2)+\frac{0.6}{2}\times(0.4x_1+0.3x_2)$

三、代码实现----Matlab

clear;clc
f = [0.4/30*2 + 0.6/2*0.4 0.4/30*5+0.6/2*0.3];
A = -ones(1,2);
b = -7;
lb = zeros(2,1);
ub = [5 6];
[x,fval] = linprog(f,A,b,[],[],lb,ub);
f1 = 2*x(1)+5*x(2);
f2 = 0.4*x(1)+0.3*x(2);

运行结果：

在这里插入图片描述

敏感性分析：

下面对结果进行敏感性分析，敏感性分析是指从定量分析的角度研究有关因素发生某种变化对某一个或一组关键指标影响程度的一种不确定分析技术。其实质是通过逐一改变相关变量数值的方法来解释关键指标受这些因素变动影响大小的规律。
改变 $f_1$ 和 $f_2$ 的权重，来观察对结果的影响。（由于两个权重和为 1 ，因此只需要改变 $f_1$ 的权重即可）

% 敏感性分析
clear;clc
A = -ones(1,2);
b = -7;
lb = zeros(2,1);
ub = [5 6];
k = (0.1:0.001:0.5);
f1 = zeros(1,length(k));
f2 = zeros(1,length(k));
y1 = zeros(1,length(k));
y2 = zeros(1,length(k));
Fval = zeros(1,length(k));
for i = 1:length(k)
    f = [k(i)/30*2 + (1-k(i))/2*0.4 k(i)/30*5+(1-k(i))/2*0.3];
    [x,fval] = linprog(f,A,b,[],[],lb,ub);
    f1(i) = 2*x(1)+5*x(2);
    f2(i) = 0.4*x(1)+0.3*x(2);
    y1(i) = x(1);
    y2(i) = x(2);
    Fval(i) = fval;
end
plot(k,f1,'g',k,f2,'-r')
legend('f1','f2')
xlabel('f1的权重')
ylabel('f1和f2的取值')
plot(k,y1,'g',k,y2,'-r')
legend('x1','x2')
xlabel('f1的权重')
ylabel('x1和x2的取值')
plot(k,Fval,'g')
legend('Fval')
xlabel('f1的权重')
ylabel('综合指标的值')

运行结果：
在这里插入图片描述

可以看出， $f_1$ 的权重的转折点在 0.333 - 0.334 之间，当 $f$ 的权重小于这个转折点时， $x_1=1,x2=6$ ，大于转折点时， $x_1=5,x_2=2$ 。其主要原因是：当 $f_1$ 的权重越小时，厂家对环境污染的权重就越大，那么厂家就更加倾向于生产污染较少的产品 B ，尽管 B 的生产费用要远高于产品 A 。

四、代码实现----python

import numpy as np
from scipy.optimize import linprog

f = np.array([0.4/30*2 + 0.6/2*0.4, 0.4/30*5+0.6/2*0.3])
A = -np.ones((1,2))
b = -7
bounds = [(0,5),(0,6)]
res = linprog(c=f,A_ub=A,b_ub=b,bounds=bounds)
x = res.x
fval = res.fun
f1 = 2*x[0]+5*x[1]
f2 = 0.4*x[0]+0.3*x[1]
print(x)
print(f1,f2,'\n',fval)

运行结果：

在这里插入图片描述

敏感性分析：

import numpy as np
from scipy.optimize import linprog
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
import mplcursors

mpl.rcParams['font.sans-serif'] = ['SimHei']   #设置字体
mpl.rcParams['axes.unicode_minus'] = False     # - 号设置

A = -np.ones((1,2))
b = -7
bounds = [(0,5),(0,6)]
k1 = 0.1
aa = []
ff1 = []
ff2 = []
X1 = []
X2 = []
Fval = []
while k1 <= 0.5:
    k2 = 1 - k1
    f = np.array([k1/30*2 + k2/2*0.4, k1/30*5+k2/2*0.3])
    res = linprog(c=f,A_ub=A,b_ub=b,bounds=bounds)
    x = res.x
    fval = res.fun
    f1 = 2.0*x[0] + 5.0*x[1]
    f2 = 0.4*x[0] + 0.3*x[1]
    aa.append(k1)    
    ff1.append(f1)
    ff2.append(f2)
    X1.append(x[0])
    X2.append(x[1])
    Fval.append(fval)
    k1 = k1 + 0.001


# 创建第一个图形
fig1, ax1 = plt.subplots()
ax1.plot(aa,ff1,'-r',label = 'f1')
ax1.plot(aa,ff2,'-g',label = 'f2')
ax1.legend()
plt.xlabel("f1的权重")
plt.ylabel("f1和f2的取值")

# 创建第二个图形
fig2, ax2 = plt.subplots()
ax2.plot(aa, X1, label='x1')
ax2.plot(aa, X2, label='x2')
ax2.legend()
plt.xlabel("f1的权重")
plt.ylabel("x1和x2的取值")

# 创建第三个图形
fig3, ax3 = plt.subplots()
ax3.plot(aa, Fval, label='fval')
ax3.legend()
plt.xlabel("f1的权重")
plt.ylabel("综合指标的值")

# 使用 mplcursors 创建交互式游标
cursor1 = mplcursors.cursor(ax1,multiple=True)
cursor2 = mplcursors.cursor(ax2,multiple=True)
cursor3 = mplcursors.cursor(ax3,multiple=True)

# 设置游标显示的内容并固定显示
@cursor1.connect("add")
def on_add(sel):
    sel.annotation.set(text=f"x={sel.target[0]:.2f}\ny={sel.target[1]:.2f}")
    sel.annotation.get_bbox_patch().set(fc="white")

@cursor2.connect("add")
def on_add(sel):
    sel.annotation.set(text=f"x={sel.target[0]:.2f}\ny={sel.target[1]:.2f}")
    sel.annotation.get_bbox_patch().set(fc="white")

@cursor3.connect("add")
def on_add(sel):
    sel.annotation.set(text=f"x={sel.target[0]:.2f}\ny={sel.target[1]:.2f}")
    sel.annotation.get_bbox_patch().set(fc="white")

plt.show()