从拉格朗日乘子--约束优化到五元齐次线性方程求解

引入
问题导入（事情起源）
三层循环
数学神器——Wolfwarm alpha
拉格朗日乘子
五元齐次线性方程
自坑三连——被伤的体无完肤
Win the final victory
番外

引入

首先我要说这是属于三个人的伟大(虽然与第三个人素未谋面，但我觉得today已经通过a math problem 产生了灵魂触碰)，在2个多小时的解题过程中，工具从百度—>CSDN—>VS17—>wolfarm alpha—>matlab，语言从数学语言—>C—>C++—>数学语言
简而言之，今晚是个奇妙之旅，下面正式起航。

问题导入（事情起源）

晚饭间，一妙人突然发来一条消息“有时间吗”，作为一直只把话言骚的基友，突然这么正经

联想最近事情，想到刚帮同学结束一场入D积极分子证明的座谈会，所以我第一反应是如下回答，但万万没想到是想我编程解决一道约束优化问题

考虑到此妙人受疫情影响电脑不在身旁，书写尚且飘逸的基础上，我又努力了一把，将以上书写体做了一下转化

求改表达式的最小值 min 0.108WA2+0.074WB2+0.055WC2+0.02WAWB+0.043WAWC+0.044WBWC0.108 W_{A}^{2}+0.074 W_{B}^{2}+0.055 W_{C}^{2}+0.02 W_{A} W_{B} +0.043 W_{A} W_{C}+0.044 W_{B} W_{C}0.108WA2+0.074WB2+0.055WC2+0.02WAWB+0.043WAWC+0.044WBWC

约束条件如下：

{WA+WB+WC=10.34WA+0.07WB+0.26WC=0.3\left\{\begin{aligned}W_{A}+W_{B}+W_{C}=1 \\ 0.34 W_{A}+0.07 W_{B}+0.26 W_{C}=0.3\end{aligned}\right.{WA+WB+WC=10.34WA+0.07WB+0.26WC=0.3

(注：公式是利用Mathpix SnipORC手写图片扫描转化）

三层循环

第一回合，直接暴力三层循环拆解。但是a,b,c精度要求是0.001，说明不是整数，此时我头脑清醒的果断投了手算一票，我预想通过两个约束条件消去目标方程2个参数，将目标方程化为一元二次方程求导，受友人“手算麻烦论”以及“电脑编程好用论”影响，更重要的是自己面子，问题始于编程，最终要是传统数学方法解决了，多没面子。

数学神器——Wolfwarm alpha

在我对代码思而无获，求而不得时，电光火石之间，我想到了传说的数学三大神器之一的Wolfwarm alpha

作为高等数学秘密武器，但我为了一道金融的题目在它身上跃跃欲试，不意外的，我撞了南墙，但我及时回了头，知道此路我已行不通。但是思路也就此停留，我迟迟无法动手开敲代码，直到两个PPT截图的到来。

拉格朗日乘子

这两张启发思维的PPT的到来，成为了今晚我们坚持不懈战胜这道题的转折点。

于是编程有了方向，灵机一动，我来抱一抱CSDN的腿，拉格朗日乘子相关内容很多，但我点开一看，不免有些失望，一来我点开几篇文章发现所讲内容重复度很高，二来都是文字理论叙述啰嗦，实战性不强，说白了，我就是想找一个实例展开并附参考代码的博文，可是没有。我想了一下可能是打开方式不对，我便直接在CSDN网页下载模块搜索拉格朗日算子线性约束，看了一下，大部分是论文级别或是课程设计级别的，并且是普遍评价都不高，浏览筛选一遍后我花15积分下载了一个资源

结果410页全是英文，并且浏览后并未发现有实例讲解，于是我开始从打开到放弃。
经历了长时间搜索我发觉这个数学题结构太简单，还没人专门编程解它，更妥当点可以说是还没人专门为他写了一篇文讲解，大概是我目前level太低，学术业务水平在淹死边缘。
就在我又一次陷入思路停留时，准备放弃时，我的友人发来了她列出的拉格朗日参数方程组，正式吹响了最后向这道题进攻的号角，此时我不得不感慨一下我的友人是一个妙人。

五元齐次线性方程

在这里，我再一次充当公式的“搬运工”

F(a,b,c,λ)=f(a,b,c)+λφ(a,b,c)+mt(a,b,c)=f(a,b,c)+λ(a+b+c−1)+m(0.34a+0.07b+0.2bc−0.3)F(a, b, c, \lambda)=f(a, b, c)+\lambda \varphi(a, b, c)+mt(a, b, c)=f(a, b, c)+\lambda(a+b+c-1)+m(0.34 a+0.07 b+0.2 b c-0.3)F(a,b,c,λ)=f(a,b,c)+λφ(a,b,c)+mt(a,b,c)=f(a,b,c)+λ(a+b+c−1)+m(0.34a+0.07b+0.2bc−0.3)

∂F∂a=0.216a+0.02b+0.043c+λ+0.34m=0\frac{\partial F}{\partial a}=0.216 a+0.02 b+0.043 c+\lambda+0.34 m=0∂a∂F=0.216a+0.02b+0.043c+λ+0.34m=0

2F∂b=0.148b+0.02a+0.0444c+λ+0.07m=0\frac{2 F}{\partial b}=0.148 b+0.02 a+0.0444 c+\lambda+0.07 m=0∂b2F=0.148b+0.02a+0.0444c+λ+0.07m=0

∂F∂r=0.11c+0.043a+0.044b+λ+0.26m=0\frac{\partial F}{\partial r}=0.11 c+0.043 a+0.044 b+\lambda+0.26 m=0∂r∂F=0.11c+0.043a+0.044b+λ+0.26m=0

∂E∂xλ=a+b+c−1=0\frac{\partial E}{\partial x^{\lambda}}=a+b+c-1=0∂xλ∂E=a+b+c−1=0

∂f∂m=0.34a+0.07b+0.2bc−0.3=0\frac{\partial f}{\partial m}=0.34 a+0.07 b+0.2 b c-0.3=0∂m∂f=0.34a+0.07b+0.2bc−0.3=0

过了一会儿，妙人说已经解出了3个未知数，此时我也有了最后编程思路，于是发起了最后进攻。

看过了这样的代码

Ar=rank(A); Br=rank(B);
B=rref(B);
if (Ar<Br)flag=0; S=0;
elseif (Ar==Br && Ar==Alie)flag=1; S=B(:,Alie+1);
else%将能构成单位矩阵的列号存储在行向量I中,即存储了极大线性无关向量的编号%将剩余列号存入行向量C中for i=1:Arfor j=1:Alieif(B(i,j)==1)I(i)=j;breakendendendC=setdiff(1:Alie,I);%由线性代数知识可得基础解系ILim=length(I); CLim=length(C);S(Alie,CLim)=0;%初始化S，S行数为A列数，S列数为C的维度for i=1:CLimS(C(i),i)=-1;for j=1:ILimS(I(j),i)=B(j,C(i));endendif(nargin==1)flag=2;elseflag=3;S(Alie,CLim+1)=0;for i=1:ArS(I(i),CLim+1)=B(i,Alie+1);endend
end

但最后代码很简单.

A=[0.216 0.02 0.043 1 0.34;0.02 0.148 0.044 1 0.07;0.043 0.044 0.11 1 0.26;1 1 1 0 0;0.34 0.07 0.26 0 0]
Z=[0 0 0 1 0.3]'
X=inv(A)*Z

但我比较不争气，在我这里又出现了出现了连续三连错误。

自坑三连——被伤的体无完肤

错误都比较弱智，在我们充斥着学术味道的聊天中，我们不断填我挖坑的过程，被完成记录下来，为读者更有身临其境之感，我将直接放图。

第一坑
(白框为挖坑者我)

在这里，我不禁害羞的低下了头，但还好，只是断了腿，我还有手和脑。

第二坑

写到这里，我觉得手已断，心已碎，只剩下脑细胞还在呼吸，于是我停下手中的工作开始做思想反思。

第三坑
在我进行深刻反思时，我的妙友找到了我的第三个也是最后一个坑，助我赢得了这场战争的胜利。

没错，我把0.26敲成0.026，结果验算结果自然失败。

此时，我不禁想要一瓶酒，一支烟，麻痹一下仅剩的神经细胞。

Win the final victory

在我纠完错误，再次运行代码得到了最终正确答案和友人手算结果相同，(没想到机算的慢一步)。在这次我要再次讲一下我的友人是一个妙人，立图为证，关于她的运算能力及准确度我要在这里夸奖一下。

今晚奇遇关于我与那位友人部分大致如上，下面放出正确结果，以示纪念意义。

番外

为故事情节的完整性，我在下面放出那位素未谋面的友人与我这位妙友之间进行的学术对话，灵魂交谈。

(注：以上部分言论是在严肃紧张中の低端学术讨论中，自我调侃语句，我们对学术是充满敬畏的)

友人说，这题是她朋友帮一位国际生写论文需要的，大概，四舍五入，我这也是为了促进国际友谊做了贡献了吧。

完