利用Matlab解决线性规划问题并绘制特定形状的空间曲面（约束区域的绘图）

今天女朋友给我发了一个问题
让我帮忙把这个空间平面y给画出来

我寻思着正好前段时间学了一些matlab的绘制曲面的方法，说不定刚好可以用上

那么就开始分析吧！

这应该就是一个高中常见的二元优化问题，但是高中碰到这种问题一般都需要一个个点去尝试，看看最后到底哪个点的值更大，这里我们有matlab绘图的话，就可以直接从图上看出来到底哪个是极值

首先我敲下了

x2=4:0.1:90;
x1=2:0.1:20;
[x1,x2]=meshgrid(x1,x2);
y=-42.9133*x1-1.9515*x2+1339.5133;
mesh(x1,x2,y)

但是不用放图大家也能想到，这其实就是一个一般的空间平面，缺少了关于x1,x2那个不等式的约束

那么该如何加上这个约束去绘图呢？

我一开始想，如果给定x1 x2的范围，

x2=4:0.1:90;
x1=2:0.1:20;

再甩上那个不等式

7.883*x1+0.862*x2>=101.314

利用matlab的solve（）函数来解一个不等式，是不是就得到了满足条件的x1 x2，然后再利用meshgrid()生成网格，就可以成功绘图了呢？

显然根据尝试的结果来看，是不行的（还临时看了看solve的用法，由于我matlab是2018版本的，solve()里的语法大变样了，我一直用老语法，总报错，在这卡了很久），根据我查的资料来看，结合尝试，是不能直接输入不等式>=去求解的，只能输入==求解方程

那么如何才可以解不等式呢？

查资料的过程中意外发现一个函数

linprog()

竟然直接就可以解决一整个线性规划的问题，我天

然后琢磨了一段时间到底咋用这个函数（具体咋用的大家可以再百度一下），我终于敲下了这几行代码，解出了女朋友发的这个问题的最优解

c=[42.9133 1.9515];
a=[-7.883 -0.862];
b=[-101.314];
lb=[2 4];
lp=[20 90];
[x,fval]=linprog(c,a,b,[],[],lb,lp);
X1=x(1)    %%%回归分析结果
X2=x(2)
Y1=-fval+1339.5133   %%(x1,x2)=(3.0138,90),y=1034.7

那行吧，现在知道最优解是啥了，还差一个图像绘制

其实也就是想要达到如下的效果：

让这个y只在约束的x1 x2范围内绘图

该咋整呢？我想了一下，这个matlab绘制三维图，是基于网格meshgrid来弄的，而且矩阵矩阵，似乎干啥都离不开矩阵这个概念，也就是说一定得基于一个给定的矩形区域才可以绘制图像，那岂不是很难弄？

于是又搜索了很多线性规划绘图啊，特定区域绘图啊，三维图像切割啥的，都没有找到满意的结果

最后终于发现一个六年前的答案

也就是说，

如果我构造一个和目标函数y同规模的全1矩阵，让深蓝色区域的矩阵元素值为0，再拿去和y逐项相乘，不就可以把y的蓝色区域清零吗！这就达到了约束区域显示的目的

Y=ones(size(y));
for i=1:861   %遍历Y的每一行，在一行中，对于前几列的元素才要清零
    for k=1:103-0.108193277*i  %设置要清零的列数的起止点
%由于我的x1 x2的间隔都是取0.1，所以可以由此计算出那条分界线的斜率和截距
        Y(i,k)=0;  
    end
end
y=y.*Y;

来看看现在的y长啥样

可以看到，约束区域以外的地方被置零了，然后利用放大镜，稍微放大一点点

就达到目的了！

然后如果想一直看到这个约束区域，可以把之前那个全1矩阵Y乘以5，再绘制出来，合理放大，就可以看到了

最后放上完整代码，欢迎交流！

clear 
clc
%%%%%%%%线性回归求解
c=[42.9133 1.9515];
a=[-7.883 -0.862];
b=[-101.314];
lb=[2 4];
lp=[20 90];
[x,fval]=linprog(c,a,b,[],[],lb,lp);
X1=x(1)    %%%输出回归分析结果
X2=x(2)
Y1=-fval+1339.5133   %%%(x1,x2)=(3.0138,90),y=1034.7
%%%%%%%%%


x1=2:0.1:20;
x2=4:0.1:90;
[x1,x2]=meshgrid(x1,x2);
y=-42.9133*x1-1.9515*x2+1339.5133;
Y=ones(size(y));  %%构造同型全1矩阵Y
for i=1:861
    for k=1:103-0.108193277*i
        Y(i,k)=0;  %%把约束区域以外的地方清零
%详细解释：Y的尺寸是861*181，对应的是x1=2:0.1:20,x1方向长18，每0.1取一个点
%所以Y有181个列，同理有861个行
%经过直线斜率和截距的计算：
%第1行的前123个元素需要清零，第953行的第1个元素需要清零
%所以可以得出，第i行的前k列元素需要清零，(953-i)/952=k/123,利用了三角形相似
    end
end
y=y.*Y;
meshz(x1,x2,y)
title('优化问题图解')
xlabel('2<=x1<=20')
ylabel('4<=x2<=90')
zlabel('y轴')
hold on
plot3(X1,X2,Y1,'g.','markersize',30) %标记处极值点
Y=Y+5; %为了加上底座，能够一直看到约束区域，把Y抬高了，这样稍微放大后不会消失
plot3(x1,x2,Y)