Yalmip入门教程(3)-约束条件的定义

        博客中所有内容均来源于自己学习过程中积累的经验以及对yalmip官方文档的翻译:https://yalmip.github.io/tutorials/

        之前的博客简单介绍了约束条件的定义方法,接下来将对其进行详细介绍。

        首先简单复习一下:

        1.定义约束条件可以使用矩阵拼接,或者’+’号,两种方式是等价的;

        2.Yalmip中的约束可以写成连续的不等式形式;

        3.Yalmip中不支持严格的不等号(<或>),一旦出现,就会报错,然后弹出“猫猫图”。

        下面将继续讲解约束条件定义的方法与技巧。

1.约束的普通形式与矩阵形式

例1:请使用Yalmip工具箱写出下列优化问题的约束条件

        采用常规的思维,可以分别写出每个约束条件,然后采用拼接的方式将所有约束条件组合起来,matlab代码如下:

x = sdpvar(3,1);C = [sum(x) <= 12 ,...5*x(1) + x(3) <= 9 ,...3*x(1) + 2*x(2) >= 6 ,...x(2) + 2*x(3) >= 1 ,...x >= 0];

        但是在求解优化问题的过程中,很多时候需要我们写出优化问题的对偶形式或KKT条件,这时候如果可以写成紧凑的矩阵形式,转换将会更加容易。例1的约束条件写成紧凑的矩阵形式为:

 

        将约束条件写作矩阵形式时,代码如下:

x = sdpvar(3,1);
A = [-1   0   0;0  -1   0;0   0  -1;1   1   1;5   0   1;-3  -2   0;0  -1  -2];
b = [0;0;0;12;9;-6;-1];
C = [A*x <= b];

        两种写法是完全等价的,下面的代码分别采用两种不同的约束定义方式对优化问题进行求解,可以得到相同的优化结果:

clc
clearx = sdpvar(3,1);
A = [-1   0   0;0  -1   0;0   0  -1;1   1   1;5   0   1;-3  -2   0;0  -1  -2];
b = [0;0;0;12;9;-6;-1];
C1 = [A*x <= b];
C2 = [sum(x) <= 12 ,...5*x(1) + x(3) <= 9 ,...3*x(1) + 2*x(2) >= 6 ,...x(2) + 2*x(3) >= 1 ,...x >= 0];
objective = [1 2 3]*x;
optimize(C1 , objective);
objective1 = value(objective)
x1 = value(x)
optimize(C2 , objective);
objective2 = value(objective)
x2 = value(x)

运行结果:

objective1 =3.3333x1 =1.33331.00000objective2 =3.3333x2 =1.33331.00000

        在实际编程过程中,两种形式各有优劣,但完全等价,可以根据自己的需要采用一般形式或矩阵形式定义约束条件,只要自己感觉更方便即可。

2.小技巧

2.1查看约束变量的参数

        在定义约束条件之后,在工作区双击对应的约束变量,或者在命令行输入约束变量再回车,可以查看约束的列表,并检查约束是否已有效定义了,例如,例1中用两种不同方式定义的约束条件列表分别如下:

        在该列表中,ID一列表示约束的编号,每个约束都是用’,’或者’+’隔开的,因此约束变量C1中仅有1个约束,约束变量C2中共有5条约束。Constraint一列表示约束的类型和数目,Element-wise inequality表示该约束条件为不等式约束,如果是等式约束,则显示‘Equality constraint’,其他一些常见的约束类型如表1所示。约束类型后面显示的是约束的数目,例如约束变量C1共包含七条约束,因此显示为7×1的约束,而约束变量C2中的第五条约束x>=0,实际上包含x1>=0,x2>=0,x3>=0共三条约束,因此显示为3×1的约束。 

表1 常见的约束类型

约束类型

含义

Matrix inequality

矩阵不等式约束

Second order cone constraint

二阶锥约束

Rotated Lorentz constraint

旋转锥约束

Binary constraint

二进制约束

Eigenvalue constraint

特征值约束

KYP constraint

Kalman Yakubovich Popov 约束

Sum-of-square constraint

平方和约束

Logic constraint

逻辑约束

Semi-continuous variable

半连续的变量约束

Semi-integer variable

半整数的变量约束

Vectorized second-order cone constraints

矢量化二阶锥约束

2.2尽量避免使用循环语句

        在Yalmip工具箱中,大部分matlab内置函数都可以用于sdpvar变量。之前我写过几篇博客说明如何在matlab代码中尽量避免循环的使用,以提升代码运行效率,因此也可将这种思想用于约束条件的定义过程中。下面以Yalmip工具箱官方文档中给出的数独求解的算例(https://yalmip.github.io/example/sudoku/)进行说明。

例2:使用Yalmip工具箱求解图1所示的数独问题

图1 一个数独问题

        为了求解这个问题,可以定义一个9×9×9的三维0-1变量A(i,j,k),当A(i,j,k)取值为1时,表示该数独的第i行第k列的取值为k,为了求解上述数独问题,需要写成初始的数独矩阵,并定义决策变量A:

A0 = [
0 4 7 0 5 0 0 0 8;
6 0 5 0 3 0 2 0 1;
0 0 0 7 0 6 0 3 0;
0 0 6 0 7 0 0 2 4;
9 0 0 8 0 4 0 0 6;
4 5 0 0 1 0 9 0 0;
0 1 0 5 0 2 0 0 0;
2 0 8 0 4 0 5 0 3;
5 0 0 0 9 0 7 1 0;
];A = binvar(9,9,9,'full');

        首先,我们需要把变量A中的元素取值和初始矩阵A0对应起来,代码如下:

for i = 1:9for j = 1:9if S(i,j)F = [F, A(i,j,A0(i,j)) == 1];endend
end

        这也是常规使用循环语句的思维,但是,如果使用find和sub2ind命令,可以避免循环语句的使用。find函数比较常用,也比较简单,这里不再介绍,重点介绍一下sub2ind。该函数用于将下标向量组转为线性索引,具体的含义我们可以通过一个例子来理解。假设在matlab中有一个矩阵3×3的矩阵S,我们如果想取出矩阵S第3行第2列的元素,可以采用下标的形式写成S(3,2),如果采用线性索引的话,也可以写成A(6)。下标和线性索引的关系如图2所示:

图2 matlab二维数组中下标和线性索引的关系

        sub2ind就可以把一组下标转为一组线性索引,以避免循环的使用,标准语法为:

ind = sub2ind(sz,row,col)

        例如,我们要取出矩阵S的元素S(1,2),S(1,3),S(2,2)和S(3,2),也就是要求出线性索引向量[4 5 6 7],采用sub2ind函数即可解决:

row = [1 2 3 1];
col = [2 2 2 3];
sz = [3 3];
ind = sub2ind(sz,row,col)

输出结果为:

ind =4     5     6     7

        因此,使用find与sub2ind,可以将上述代码改写成:

[i,j,k] = find(A0);
F = [F, A(sub2ind(size(A),i,j,k)) == 1];

        这样做就可以避免使用双层循环。另外,数独游戏的规则是,每一行、每一列以及每一个3×3的子块的数字不能重复,其中每行每列数字不能重复的约束比较简单,代码如下:

F = [sum(A,1) == 1, sum(A,2) == 1, sum(A,3) == 1];

        要保证每个3×3的子块数字不能重复,可以采用循环的方式,代码如下:

for m = 1:3for n = 1:3for k = 1:9s = sum(sum(A((m-1)*3+(1:3),(n-1)*3+(1:3),k)));             F = [F, s == 1];endend
end

        为了尽量避免使用循环语句,可以将代码改写如下:

for k = 1:9F = [F, sum(sum(A(i1:i1+2 , j1:j1+2 , :))) == 1];
end

        这样写可以只用一层循环,具体原理比较简单,这里不再过多讲解。另外,对于数独问题,我们只需要求出可行解,并没有优化目标,因此使用optimize函数时可以只输入约束条件。

sol = optimize(F);

        求解成功后,可以使用一个矩阵Z来存储最终的数独求解结果:

Z = 0;
for i = 1:9Z = Z  + i*value(A(:,:,i));
end
Z

        为了避免循环语句,也可以写成:

Z = value(reshape(A,9,81)*kron((1:9)',eye(9)))

        其中,reshape函数用于改变矩阵的形状,并返回一个新的矩阵。在这条语句中,reshape(A,9,81)表示将矩阵A重塑为一个9行81列的矩阵。kron函数用于计算两个矩阵的张量积,在这条语句中,kron((1:9)',eye(9))可以得到一个81行9列的矩阵,其中对应位置的元素由(1:9)'和单位矩阵的元素相乘得到。将矩阵reshape(A,9,81)与矩阵kron((1:9)',eye(9)),得到的就是9×9的矩阵,也就是数独问题的解。

        综上所示,使用常规思维+更多的循环语句求解数独问题的完整代码如下:

方法1:放肆使用循环语句

clc
cleartic
A0 = [
0 4 7 0 5 0 0 0 8;
6 0 5 0 3 0 2 0 1;
0 0 0 7 0 6 0 3 0;
0 0 6 0 7 0 0 2 4;
9 0 0 8 0 4 0 0 6;
4 5 0 0 1 0 9 0 0;
0 1 0 5 0 2 0 0 0;
2 0 8 0 4 0 5 0 3;
5 0 0 0 9 0 7 1 0;
];A = binvar(9,9,9,'full');
F = [sum(A,1) == 1, sum(A,2) == 1, sum(A,3) == 1];for i = 1:9for j = 1:9if A0(i,j)F = [F, A(i,j,A0(i,j)) == 1];endend
endfor m = 1:3for n = 1:3for k = 1:9s = sum(sum(A((m-1)*3+(1:3),(n-1)*3+(1:3),k)));             F = [F, s == 1];endend
enddiagnostics = optimize(F);Z = 0;
for i = 1:9Z = Z  + i*value(A(:,:,i));
end
Z
toc

方法2:尽量减少循环语句

clc
cleartic
A0 = [
5 3 0 0 7 0 0 0 0;
6 0 0 1 9 5 0 0 0;
0 9 8 0 0 0 0 6 0;
8 0 0 0 6 0 0 0 3;
4 0 0 8 0 3 0 0 1;
7 0 0 0 2 0 0 0 6;
0 6 0 0 0 0 2 8 0;
0 0 0 4 1 9 0 0 5;
0 0 0 0 8 0 0 7 9;
];A = binvar(9,9,9,'full');F = [sum(A,1) == 1, sum(A,2) == 1, sum(A,3) == 1];i1 = [1,1,1,4,4,4,7,7,7];
j1 = [1,4,7,1,4,7,1,4,7];
for k = 1:9F = [F, sum(sum(A(i1:i1+2 , j1:j1+2 , :))) == 1];
end[i,j,k] = find(A0);
F = [F, A(sub2ind(size(A),i,j,k)) == 1];diagnostics = optimize(F);Z = value(reshape(A,9,81)*kron((1:9)',eye(9)))
toc

        两种方式的运行结果完全一致:

Z =3     4     7     2     5     1     6     9     86     8     5     4     3     9     2     7     11     2     9     7     8     6     4     3     58     3     6     9     7     5     1     2     49     7     1     8     2     4     3     5     64     5     2     6     1     3     9     8     77     1     3     5     6     2     8     4     92     9     8     1     4     7     5     6     35     6     4     3     9     8     7     1     2

        但在运行时间上,两种方法存在差异,某次测试中,方法1的运行时间为0.638127秒,方法2的运行时间为0.473180秒,避免循环的使用提高了25.85%的运行效率。

        减少循环语句的使用,在小规模问题中效率提升可能不是很明显,但如果优化问题比较复杂,约束也很复杂,如果可以把尽量减少循环语句的思想应用到实际编程中,效率提升可能是非常大的。虽然思考如何减少循环语句可能要废很多脑细胞,但我觉得应该会挺值得。(我自己经历过一个优化问题求解需要8个小时左右,每次等待结果、调试代码都非常痛苦,将近半个月都没有调好,后来想通了,花了两三天的时间优化代码,把运行时间缩短到了1个半小时左右,很快就调试成功了)

3.测试题

3.1测试1

使用矩阵形式的约束条件求解下列优化问题:

max z=3x1+x2

s.t. x1-x2≥-2

x1-2x2≤3

3x1+2x2≤14

x1≥0, x2≥0

3.2测试2

使用Yalmip工具箱求下列数独问题的解:

3.3测试3

        改写下面的代码以避免使用循环语句,并对比两种方式的代码运行时间。

clc
clear
yalmip('clear')tic
n = 50;
a = rand(n,n,n);
x = sdpvar(n,n,n);
z = 0;
for k = 1:nfor kk = 1:nfor kkk = 1:nz = z + a(k,kk,kkk) * x(k,kk,kkk);endend
end
toc

3.4测试题参考答案

        第三章测试题的参考答案可以从下面的链接中获取:

https://download.csdn.net/download/weixin_44209907/88062937

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/4301.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

GRE和MGRE

目录 GRE GRE环境的搭建 MGRE MGRE的配置 MGRE环境下的RIP网络 MGRE实验 VPN 说到GRE&#xff0c;我们先来说个大家熟悉一点的&#xff0c;那就是VPN技术。 背景需求 企业、组织、商家等对专用网有强大的需求。 高性能、高速度和高安全性是专用网明显的优势。 物理专…

Notepad++ 配置python虚拟环境(Anaconda)

Notepad配置python运行环境步骤&#xff1a; 打开Notepad ->”运行”菜单->”运行”按钮在弹出的窗口内输入以下命令&#xff1a; 我的conda中存在虚拟环境 (1) base (2) pytorch_gpu 添加base环境至Notepad中 cmd /k chdir /d $(CURRENT_DIRECTORY) & call cond…

TX Barcode .NET for WPF Crack

TX Barcode .NET for WPF Crack 用于WPF软件的TX Barcode.NET包括一天完成的功能以及用于WPF的软件的2D条形码控制。 用于WPF的TX Barcode.NET的功能和属性&#xff1a; 它具有以下特性和属性&#xff0c;如&#xff1a; 常见的文字处理功能&#xff1a;它可以为用户和开发人员…

FastEdit ⚡:在10秒内编辑大型语言模型

概述&#xff1a; 这个仓库旨在通过一个单一的命令&#xff0c;有效地将新鲜且定制化的知识注入到大型语言模型中&#xff0c;以辅助开发人员的工作。 支持的模型&#xff1a;○ GPT-J (6B)○ LLaMA (7B/13B)○ BLOOM (7.1B)○ Falcon (7B)○ Baichuan (7B/13B)○ InternLM (7…

stable diffusion windows本地搭建的坑

刚刚2小时前&#xff0c;我搭好了&#xff0c;欣喜若狂&#xff0c;开放端口&#xff0c;同事也尝试了。我的配置 16G内存&#xff0c;AMD卡&#xff0c;有gpu但是没有用。这里不说具体步骤&#xff0c;只说坑点。 首先就是安装gfpgan、clip、openclip问题&#xff0c;我参考…

Kind | Kubernetes in Docker 把k8s装进docker!

有点像杰克船长的黑珍珠 目录 零、说明 一、安装 安装 Docker 安装 kubectl 安装 kind 二、创建/切换/删除集群 创建 切换 删除 将镜像加载到 kind 群集中 零、说明 官网&#xff1a;kind Kind&#xff1a; Kubernetes in Docker 的简称。kind 是一个使用 Docker 容…

CentOS5678 repo源 地址 阿里云开源镜像站

CentOS5678 repo 地址 阿里云开源镜像站 https://mirrors.aliyun.com/repo/ CentOS-5.repo https://mirrors.aliyun.com/repo/Centos-5.repo [base] nameCentOS-$releasever - Base - mirrors.aliyun.com failovermethodpriority baseurlhttp://mirrors.aliyun.com/centos/$r…

数据仓库建设-数仓分层

数据仓库能够帮助企业做出更好的决策&#xff0c;提高业务效率和效益&#xff1b;在数据仓库建设时&#xff0c;绕不开的话题就是数仓分层。 一、数据分层的好处 1. 降低数据开发成本 通用的业务逻辑加工好&#xff0c;后续的开发任务可以基于模型快速使用&#xff0c;数据需…

IDEA使用方式

1.翻译 1.Plugins插件&#xff1a;Chinese中文插件 文件F 编辑E 视图V 导航N 代码C 分析Z 重构R 构建B 运行U 工具T VCSS 窗口W 帮助H文件N 新建N 打开 打开最近 关闭项目 设置T 项目结构 文件属性 保存全部S 从磁盘全部重新加载 作废缓存/重启 导出/导入操作 其他设置 导出 打…

windows下mysql8定时备份,bat脚本编写,dos免密执行

前提&#xff1a;mysql8已经安装。 编写脚本copy_mysql_data.bat echo off set timestamp%date:~0,4%%date:~5,2%%date:~8,2%_%time:~0,2%%time:~3,2%%time:~6,2% set backupfileD:\ProgramData\MySQL\Backup\backup_%timestamp%.sql set mysqlpathD:\Program Files\MySQL\MyS…

ubuntu版本Linux操作系统上安装键盘中文输入法

要在ubuntu版本Linux操作系统上安装键盘中文输入法 可以按照以下步骤进行操作&#xff1a; 1、Linux终端输入&#xff1a;sudo apt-get install ibus-pinyin 这将安装一个常用的中文输入法 “ibus-pinyin”。 2、重新启动系统&#xff1a;为了使输入法生效&#xff0c;需要…

ubuntu打开usb摄像头

文章目录 前言一、识别 usb 摄像头二、安装应用程序显示摄像头捕捉到的视频1、使用应用程序茄子&#xff08;cheese&#xff09;2、运行 cheese 捕捉视频 前言 记录一下解决在 Linux 下打开 usb 摄像头界面黑屏的问题。 一、识别 usb 摄像头 1、保持在 ubuntu 界面&#xff0…

前端学习记录~2023.7.17~CSS杂记 Day9

前言一、浮动1、使盒子浮动起来2、清除浮动3、清除浮动元素周围的盒子&#xff08;1&#xff09;clearfix 小技巧&#xff08;2&#xff09;使用 overflow&#xff08;3&#xff09;display: flow-root 二、定位1、定位有哪些2、top、bottom、left 和 right3、定位上下文4、介绍…

科技赋能企业,实现数字化转型

科技是第一生产力&#xff0c;数字技术即科技&#xff0c;可以改变传统的商业模式&#xff0c;为各行各业注入新的活力。 推动企业数字化转型&#xff0c;可是实现行业的效率提升&#xff0c;实现跨界重组&#xff0c;重构产业模式&#xff0c;为产业格局重新赋能&#xff0c;最…

AJAX: 事件循环(举例细论)

概念&#xff1a;执行任务和收集异步任务&#xff0c;在调用栈空闲时&#xff0c;反复调用任务队列里回调函数的一种执行机制 原因&#xff1a;JavaScript 是单线程的&#xff0c;为了不阻塞 JS 引擎&#xff0c;设计执行代码的模型 JS内代码如何执行&#xff1a; 执行同步代…

PHP与Golang对战:两种语言的比较与应用场景探讨

引言 在软件开发领域&#xff0c;选择一种合适的编程语言对于项目的成功至关重要。而在今天的文中&#xff0c;我们将探讨两个备受争议的编程语言——PHP与Golang之间的对战。通过比较它们的优势和应用场景&#xff0c;帮助开发者更好地了解如何选择适合自己项目的语言。 PHP的…

伙伴云CEO戴志康:我们为什么要做伙伴云?

分享嘉宾&#xff1a;戴志康&#xff0c;伙伴云CEO 以下为演讲实录⬇⬇⬇ 01选择人更少的一条路&#xff0c;从B级走向A级 我一直想和大家交流一个话题&#xff0c;关于我们为什么要做伙伴云。既代表我自己&#xff0c;同时也代表我们团队的一些想法。 我是一个怀疑论者。大…

Win10 配置NDK安装2023.7.19版本

NDK安装流程 1. 下载&#xff1a;2. 安装&#xff1a;3. 测试&#xff1a; 在大多数情况下&#xff0c;使用 Android SDK 管理器安装 NDK 会更轻松。本文单独安装NDK&#xff0c;但后续也可以使用管理器进行管理。 1. 下载&#xff1a; 地址 Fig.1 最新稳定版本 2. 安装&…

3.6 Bootstrap 导航元素

文章目录 Bootstrap 导航元素表格导航或标签胶囊式的导航菜单基本的胶囊式导航菜单垂直的胶囊式导航菜单 两端对齐的导航禁用链接下拉菜单带有下拉菜单的标签带有下拉菜单的胶囊标签页与胶囊式标签页 Bootstrap 导航元素 本文将讲解 Bootstrap 提供的用于定义导航元素的一些选项…

OpenCv之图像直方图

目录 一、基本概念 二、使用OpenCv统计直方图 三、使用掩膜的直方图 一、基本概念 图像直方图是用一表示教字图像中亮度分布的直方图&#xff0c;标绘了图像中每个高度值的像素数。可以借助观察该有方图了解需要如何调整亮度分布的直方图。这种直方图中&#xff0c;横坐标的左…