Matlab求解优化问题-Yalmip

1.前言

       前述文章介绍fmincon函数用以求解非线性规划问题，但是fmincon面对非线性整数规划没有办法。对此，根据对网上资源的搜索，发现有一个工具箱可以解决这一个问题，那就是Yalmip。

       Yalmip是一个建模工具，也可以称为一种“语言”，通过这种“语言”来描述模型，然后再调用其他求解器（如gurobi、cplex等）来求解模型。相当于一个将“yalmip语言”转换成其他求解器“语言”的语言转换器。

       Yalmip的优势：

• Yalmip可以使用各种求解器，包括clpex,scip等
• 在Yalmip中决策变量可以传入矩阵

1.1 Yalmip安装

• 官网下载https://yalmip.github.io/
• 解压至matlab/toolbox（自己记住Matlab安装到哪哦！）
• matlab中设置路径
  • 主页-设置路径-找到Yalmip文件夹-添加
• 检查是否配置成功：matlab中调用yalmitest命令，查看所有支持的求解器已经他们的安装状态。出现以下界面说明配置成功。

1.2 Yalmip使用

• 变量
  • sdqvar()设置实型变量，
  • intvar()设置整形变量，
  • binvar()设置0-1变量。
  • 示例：x=binvar(20,92,‘full’);%创建0/1型决策变量;full表示全参数矩阵
• 约束
  • 可以先设置一个空矩阵，然后往里面添加约束就好了
  • 示例：

C=[x(15,28)==1,x(15,29)==1,x(16,38)==1,x(2,39)==1,x(7,61)==1,x(20,92)==1];%特殊点分配，6个约束
%每个路口节点j被分配到1个平台i，92个约束
for j=1:92
    C=[C,sum(x,1)==1];%sum函数2对每列的数相加，输出行(j)
end

• 目标函数
  • 示例：$Z=max(x\ast c{n}_j)-min(x\ast c{n}_j)$;
• 求解
  • 使用sdpsettings()函数设置参数，这个函数里面实际上可以写参数也可以不写参数
  • ops=sdpsettings();
  • 求解使用optimize()函数，其中optimize(C,Z,ops)中的C是约束条件，Z是目标函数，ops是参数设置。
  • 使用double()函数，将求解的实型变量x1、x2以及目标函数值f转换为实数。
  • z=double(Z);x=double(x);
• 官方文档：https://yalmip.github.io/tutorial/basics/

2.模型和代码

2.1 数据

       我们仍然以2011年高教社杯全国大学生数学建模竞赛Ｂ题—交警服务平台的设置与调度提供的数据为例进行建模与分析。

       数据介绍：该数据集给出了某个城市A城区的交通路口以及交巡警平台信息，包括路口位置、路线的起点和终点路口节点、路口节点的发案率（每个路口平均每天的发生报警案件数量）、交巡警平台所在路口节点以及出入各个区的路口节点。该城区有92个路口节点（记为1-92），20个交巡警路口平台（记为1-20）。

2.2 假设、符号

       此例子涉及的假设有：

• 交巡警警车在市区道路的行驶速度恒为60km/h（虽然现实并非如此）；
• 在报警事件发生时，交巡警应尽量在3分钟之内到达。

表2 符号定义表

符号	含义
m	交巡警服务平台个数，例子中m= 20
n	交通路口节点个数，例子中n= 92
v	交巡警警车在市区道路的行驶速度，v=60 km/h
$c{n}_j$	第j个节点事件发生率（每个路口平均每天的发生报警案件数量）
$x_{ij}$	第i个平台是否管辖第j个节点,0-1变量
$w_{ij}$	最短距离是否小于最长响应时间对应的距离，例子中为3km

注：已经使用python的networkx包all_pairs_dijkstra_path_length函数计算交巡警服务平台到各个路口节点的最短距离（此函数背后的算法是Dijkstra算法）。由于篇幅限制，有需要的朋友可发送“shortest_path_distance’到后台获取代码。

2.3 模型

       模型原则：尽量满足3分钟到达报警现场，对于一些特殊节点（所有平台都无法做到3分钟到达，也就是所有平台到该节点的最短距离都大于3km），安排离它们最近的服务平台管辖,以便尽快响应报警事件。

表2 特殊路口节点管辖安排

交通路口节点标号	归属的服务平台标号	最近距离（km）
28	15	4.75
29	15	5.70
38	16	3.41
39	2	3.68
61	7	4.19
92	20	3.60

       接下来，我们以最大工作量与最小工作量之差最小化为目标建立数学模型，具体模型如下：
       $minz=max(q_i)-min(q_i)(1)$
       $s.t.q_i={\sum }_{j=1}^n{x}_{ij}c{n}_j,i\in I_1(2)$
       $x_{15,28}=x_{15,29}=x_{16,38}=x_{2,39}=x_{7,61}=x_{20,92}=1(3)$
       ${\sum }_{i=1}^m{x}_{ij}=1,j\in J_1(4)$
       ${\sum }_{i=1}^m{w}_{ij}{x}_{ij}=1,j\notin (28,29,38,39,61,92)(5)$
       $x_{ij}=0或1,i\in I_1,j\in J_1(6)$
       $I_1=1,2,...,20(7)$
       $J_1=1,2,...,92(8)$

      约束（2）计算i平台的工作量，式子（3）记录特殊节点的分配结果，式子（4）保证每个路口j都有其管辖的平台，式子（5）保证非特殊交通路口节点能被其管辖平台3分钟到达，式子（6）说明$x_{ij}$为0-1变量。

2.4 代码

clear
%1.读取数据
%读取wij数据
opts = detectImportOptions('C:\Users\xxx\Desktop\distance.xls');%使用一个变量创建 SpreadsheetImportOptions 对象
opts.Sheet = 'sheet3';%读取的sheet_name
opts.SelectedVariableNames=(2:93);%读取指定列，此处为读取第2列到93列
opts.DataRange= '2:21';%读取行
w=readmatrix('C:\Users\xxx\Desktop\distance.xls',opts);%读取w矩阵
%读取cn_j数据
opts1 = detectImportOptions('C:\Users\xxx\Desktop\crimedata.xls');%使用一个变量创建 SpreadsheetImportOptions 对象
opts1.Sheet = 'sheet1';%读取的sheet_name
opts1.SelectedVariableNames=[2];%读取指定列，此处为读取第2列
opts1.DataRange= '2:93';
cn_j=readmatrix('C:\Users\xxx\Desktop\crimedata.xls',opts1);%读取cn矩阵

%2.建立模型
%2.1声明变量
x=binvar(20,92,'full');%创建0/1型决策变量;full表示全参数矩阵

%2.2目标函数
Z=max(x*cn_j)-min(x*cn_j);

%2.3约束条件
C=[x(15,28)==1,x(15,29)==1,x(16,38)==1,x(2,39)==1,x(7,61)==1,x(20,92)==1];%特殊点分配

%每个路口节点j被分配到1个平台i
for j=1:92
    C=[C,sum(x,1)==1];%sum函数2对每列的数相加，输出行(j)
end
%覆盖约束
for j=1:92
    C=[C,sum(w.*x,1)==1];%j不等于特殊节点
end

%3.参数设置和求解
ops = sdpsettings('verbose',2,'debug',1); %语法规则其实是‘field',value，不指定求解器，verbose展示求解细节的设置。0表示完全不显示，1表示适度显示，2则是完全显示；debug：当设置为1时，yalmip会将出错的原因和位置显示在命令行窗口
sol  = optimize(C,Z,ops);% sol = solvesdp(Constraints,Z,ops);
if sol.problem== 0
    xp=value(x);
    Zp=value(Z);
    qi=round(xp*cn_j,2);
else
    disp('求解过程中出错');
end