前言
方便科研人员自助建立拉格朗日模型，并直接生成对应矩阵，方便仿真（matlab版本高于2020b）。
本文基于我的上一篇文章进行了集成工作，修复了一些BUG。现在可以直接建立多自由度的模型，理论上对自由度的数量没有限制。
一、操作演示
现在我想通过拉格朗日方法建立一个数学模型，模型的状态量维度为5，并计算好了模型对应的拉格朗日算子L。
调用如下函数：
num=5;
[matrix_M,matrix_C,vector_G] = modle_builder(L,num)
运行结果：
matrix_M =
5×5 string 数组
"(p1*p3^2)/3 + (p2*p3^2)/4" "(p2*p3^2*cos(n1(t) + n2(t)))/4" "-(p2*p3*p4*sin(n1(t) - n3(t)))/4" "(p2*p3*p4*sin(n1(t) + n4(t)))/4" "-(p2*p3*p6*cos(n1(t) - n5(t)))/2"
"(p2*p3^2*cos(n1(t) + n2(t)))/4" "(p1*p3^2)/3 + (p2*p3^2)/4" "(p2*p3*p4*sin(n2(t) + n3(t)))/4" "-(p2*p3*p4*sin(n2(t) - n4(t)))/4" "-(p2*p3*p6*cos(n2(t) + n5(t)))/2"
"-(p2*p3*p4*sin(n1(t) - n3(t)))/4" "(p2*p3*p4*sin(n2(t) + n3(t)))/4" "(p2*p4^2)/4" "-(p2*p4^2*cos(n3(t) + n4(t)))/4" "-(p2*p4*p6*sin(n3(t) - n5(t)))/2"
"(p2*p3*p4*sin(n1(t) + n4(t)))/4" "-(p2*p3*p4*sin(n2(t) - n4(t)))/4" "-(p2*p4^2*cos(n3(t) + n4(t)))/4" "(p2*p4^2)/4" "-(p2*p4*p6*sin(n4(t) + n5(t)))/2"
"-(p2*p3*p6*cos(n1(t) - n5(t)))/2" "-(p2*p3*p6*cos(n2(t) + n5(t)))/2" "-(p2*p4*p6*sin(n3(t) - n5(t)))/2" "-(p2*p4*p6*sin(n4(t) + n5(t)))/2" "p2*p6^2"
matrix_C =
5×5 string 数组
"0" "(- (p2*p3^2*cos(n1(t))*sin(n2(t)…" "((p2*p3*p4*cos(n1(t))*cos(n3(t))…" "((p2*p3*p4*cos(n1(t))*cos(n4(t))…" "((p2*p3*p6*cos(n1(t))*sin(n5(t))…"
"(- (p2*p3^2*cos(n1(t))*sin(n2(t)…" "0" "((p2*p3*p4*cos(n2(t))*cos(n3(t))…" "((p2*p3*p4*cos(n2(t))*cos(n4(t))…" "((p2*p3*p6*cos(n2(t))*sin(n5(t))…"
"(- (p2*p3*p4*cos(n1(t))*cos(n3(t…" "((p2*p3*p4*cos(n2(t))*cos(n3(t))…" "0" "((p2*p4^2*cos(n3(t))*sin(n4(t)))…" "((p2*p4*p6*cos(n3(t))*cos(n5(t))…"
"((p2*p3*p4*cos(n1(t))*cos(n4(t))…" "(- (p2*p3*p4*cos(n2(t))*cos(n4(t…" "((p2*p4^2*cos(n3(t))*sin(n4(t)))…" "0" "((p2*p4*p6*sin(n4(t))*sin(n5(t))…"
"((p2*p3*p6*cos(n5(t))*sin(n1(t))…" "((p2*p3*p6*cos(n2(t))*sin(n5(t))…" "(- (p2*p4*p6*cos(n3(t))*cos(n5(t…" "((p2*p4*p6*sin(n4(t))*sin(n5(t))…" "0"
vector_G =
5×1 string 数组
"(p3*p8*cos(n1(t))*(p1 + p2))/2"
"(p3*p8*cos(n2(t))*(p1 + p2))/2"
"(p2*p4*p8*sin(n3(t)))/2"
"(p2*p4*p8*sin(n4(t)))/2"
"-p2*p6*p8*cos(n5(t))"
>>
这里矩阵M、C和向量G的含义就是最标准的矩阵方程，如下：
M q ¨ + C q ˙ + G = F M\ddot q + C\dot q + G = F Mq¨+Cq˙+G=F
二、使用步骤
matlab文件中包含多个程序，其中一个主程序和多个辅助程序，下面分别介绍
1.Matlab主程序（modle_builder）
代码如下：
function [matrix_M,matrix_C,vector_G] = modle_builder(L,num)
syms n1(t) n2(t) n3(t) n4(t) n5(t) n6(t) n7(t) n8(t) n9(t) n10(t)
syms dn1(t) dn2(t) dn3(t) dn4(t) dn5(t) dn6(t) dn7(t) dn8(t) dn9(t) dn10(t) dn11(t)
syms ddn1(t) ddn2(t) ddn3(t) ddn4(t) ddn5(t) ddn6(t) ddn7(t) ddn8(t) ddn9(t) ddn10(t)
syms p1 p2 p3 p4 p5 p6 p7 p8 p9 h1 h2 h3 h4 h5 h6 h7 h8 h9
dn=["dn1(t)","dn2(t)","dn3(t)","dn4(t)","dn5(t)","dn6(t)","dn7(t)","dn8(t)","dn9(t)","dn10(t)"];
n=["n1(t)","n2(t)","n3(t)","n4(t)","n5(t)","n6(t)","n7(t)","n8(t)","n9(t)","n10(t)"];
eq=["a"];
for i=1:num
eq(i)=string(simplify(str2sym(replace_diff(char(diff(diff(L,str2sym(dn(i))),t))))-diff(L,str2sym(n(i)))));
end
[matrix_M,matrix_C,vector_G]=element(eq);
end
2.Matlab辅助程序1（element）
主程序中调用了多个辅助程序，其中最关键的就是element.
代码如下（示例）：
function [matrix_M,matrix_C,vector_G] = element(eq)
obj =match.Yang_V1;
%%%%%%%%%%%%%%%需要根据模型修改的部分%%%%%%%%%%%%%%%
syms n1(t) n2(t) n3(t) n4(t) n5(t) n6(t) n7(t) n8(t) n9(t) n10(t)
syms dn1(t) dn2(t) dn3(t) dn4(t) dn5(t) dn6(t) dn7(t) dn8(t) dn9(t) dn10(t) dn11(t)
syms ddn1(t) ddn2(t) ddn3(t) ddn4(t) ddn5(t) ddn6(t) ddn7(t) ddn8(t) ddn9(t) ddn10(t)
syms p1 p2 p3 p4 p5 p6 p7 p8 p9 h1 h2 h3 h4 h5 h6 h7 h8 h9
ddn=["ddn1(t)","ddn2(t)","ddn3(t)","ddn4(t)","ddn5(t)","ddn6(t)","ddn7(t)","ddn8(t)","ddn9(t)","ddn10(t)"];
dn=["dn1(t)","dn2(t)","dn3(t)","dn4(t)","dn5(t)","dn6(t)","dn7(t)","dn8(t)","dn9(t)","dn10(t)"];
dnf=["dn1(t)^2","dn2(t)^2","dn3(t)^2","dn4(t)^2","dn5(t)^2","dn6(t)^2","dn7(t)^2","dn8(t)^2","dn9(t)^2","dn10(t)^2"];
%%%%%%%%%%%%%%%%定义容器存储运动方程%%%%%%%%%%%%%%%%
temp_eq=eq;
%%%%%%%%%%%%%%%%%%%%%%辅助变量%%%%%%%%%%%%%%%%%%%%%
i_lastf=size(temp_eq);
i_lastl=i_lastf(2);
j_lastf=size(eq);
j_lastl=j_lastf(2);
%%%%%%%%%%%%初始化矩阵,类型确定为字符串%%%%%%%%%%%%%
matrix_M=["a"];
matrix_C=["a"];
%%%%%%%%%%%%%%%%%%%%矩阵M赋值%%%%%%%%%%%%%%%%%%%%%%
for i = 1:i_lastl %%外循环，方程
for j =1:j_lastl %%内循环，变量
coef=coeffs(str2sym(temp_eq(i)),str2sym(ddn(j))); %%直接求系数
siz=size(coef); %%获取系数向量维度
if(siz(2)~=1) %%系数向量维度不等于1，则一定存在二阶变量对应的系数，且一定位于coef(2)
matrix_M(i,j)=string(coef(2)); %%直接赋值
temp_eq(i)=string(simplify(str2sym(temp_eq(i))-coef(2)*str2sym(ddn(j))));%%减去当前方程中对应的二阶项
elseif((siz(2)==1)&&(contains(temp_eq(i),ddn(j)))==1) %%如果方程此时仅包含该二阶变量，则此时的向量维度等于1，考虑到这种情况（一般来说不会存在）
matrix_M(i,j)=string(coef(1));
temp_eq(i)=string(simplify(str2sym(temp_eq(i))-coef(1)*str2sym(ddn(j))));
else
matrix_M(i,j)=0; %%其他情况，则没有该二阶变量。
end
end
end
%%此时temp_eq中没有了二阶状态量对应的项（除了扰动）%%
for i = 1:i_lastl %%外循环
for j =1:j_lastl %%内循环
temp1=coeffs(expand(str2sym(temp_eq(i))),str2sym(dn(j))); %%用temp1储存第i个方程对第j个速度量求的系数向量
siz=size(temp1); %%记录系数向量的维度
if((siz(2)==1)&&(contains(temp_eq(i),dn(j)))~=1) %%维度是1，相关项
matrix_C(i,j)=0;
elseif((siz(2)==1)&&(contains(temp_eq(i),dnf(j)))==1) %%维度是1，平方项的系数
matrix_C(i,j)=['(',char(temp1(1)),')','*',char(dn(j))];
elseif((siz(2)==1)&&(contains(temp_eq(i),dnf(j)))~=1&&(contains(temp_eq(i),dn(j)))==1) %%维度是1，但是是一次项的系数
matrix_C(i,j)=['(',char(temp1(1)),')','*',char(dn(j))];
elseif((siz(2)==2)&&(contains(temp_eq(i),dnf(j)))~=1) %%维度是2，不包含平方项，则第二个是一次项的系数，直接处理赋值
matrix_C(i,j)=char(temp1(2));
matrix_C(i,j)=string(expand(str2sym(string(matrix_C(i,j)))));
matrix_C(i,j)=replace_point(matrix_C(i,j));
matrix_C(i,j)=obj.getValue(matrix_C(i,j),j,j_lastf(2));
matrix_C(i,j)=string(simplify(str2sym(string(matrix_C(i,j)))));
elseif((siz(2)==2)&&(contains(eq(i),dnf(j)))==1) %%维度是2，包含平方项，，这里假设了这种情况，肯定是有常数的，即常数+平方的组合，没考虑一次和平方的组合
matrix_C(i,j)=['(',char(temp1(2)),')','*',char(dn(j))];
elseif(siz(2)==3) %%包含一次项和平方项，则分别处理
matrix_C(i,j)=char(temp1(2));
matrix_C(i,j)=replace_point(matrix_C(i,j));
matrix_C(i,j)=obj.getValue(matrix_C(i,j),j,j_lastf(2));
matrix_C(i,j)=string(expand(str2sym(char([char(matrix_C(i,j)),'+',' ','(',char(temp1(3)),')','*',char(dn(j))]))));
matrix_C(i,j)=string(simplify(str2sym(string(matrix_C(i,j)))));
end
end
end
matrix_C_temp=str2sym(matrix_C);
for i=1:i_lastl
for j=1:j_lastl
temp_eq(i)=string(simplify(str2sym(temp_eq(i))-matrix_C_temp(i,j)*str2sym(dn(j))));
end
end
vector_G=temp_eq';
end
3.Matlab其他辅助程序
这里包含一些字符串的替换程序
代码如下：
function [y1] = replace_t(x1)
y1=replace(x1, '(t)', '');
end
function [y1] = replace_quad(x1)
y1=char(replace(x1, ' ', ''));
end
function [y1] = replace_point(x1)
replace(x1,"1.0*","");
y1=replace(ans,"2.0*","2*");
end
function [y1] = replace_diff(x1)
replace(x1, 'diff(n1(t), t)', 'dn1(t)');
replace(ans, 'diff(n2(t), t)', 'dn2(t)');
replace(ans, 'diff(n3(t), t)', 'dn3(t)');
replace(ans, 'diff(n4(t), t)', 'dn4(t)');
replace(ans, 'diff(n5(t), t)', 'dn5(t)');
replace(ans, 'diff(n6(t), t)', 'dn6(t)');
replace(ans, 'diff(dn1(t), t)', 'ddn1(t)');
replace(ans, 'diff(dn2(t), t)', 'ddn2(t)');
replace(ans, 'diff(dn3(t), t)', 'ddn3(t)');
replace(ans, 'diff(dn4(t), t)', 'ddn4(t)');
replace(ans, 'diff(dn5(t), t)', 'ddn5(t)');
y1=replace(ans, 'diff(dn6(t), t)', 'ddn6(t)');
end
function [y1] = replace_point(x1)
replace(x1,"1.0*","");
y1=replace(ans,"2.0*","2*");
end
4.最后一步，JAR包调用
因为文件中调用了jar包中的一个静态方法，这里直接的静态方法如下：
这里不熟悉的人可能有点懵，但实际上就是编写了一段处理字符串的JAVA方法，直接建立好对应的包，复制程序，导出jar包，最后将jar包粘贴至MATLAB对应的位置。
还是不理解的可以自行搜索“MATLAB如何调用jar包”
package match;
public class Yang_V1 {
public static void main(String[] args) {
}
public static String getValue(String str,int ingex,int number){
String[] s = str.split(" ");
String outa="";
if(s.length%2==1){
String[] temp = new String[(s.length+1)/2];
for (int i = 0; i < temp.length; i++) {
if(i==0){
temp[0]=s[0];
}
else{
temp[i]=s[i*2-1]+s[2*i];
}
}
for (int i = 0; i < temp.length; i++) {
for (int i1 = 0; i1 < number; i1++) {
if((temp[i].contains("dn"+(i1+1)+"(t)")&&((i1+1)!=ingex))){
if (temp[i].startsWith("-2")){
temp[i]="-"+temp[i].substring(3,temp[i].length());
} else if (temp[i].startsWith("2")) {
temp[i]=temp[i].substring(2,temp[i].length());
}
}
}
}
for (int i = 0; i < temp.length; i++) {
outa =outa+temp[i];
}
}else {
String[] temp = new String[(s.length)/2];
for (int i = 0; i < temp.length; i++) {
temp[i]=s[i*2]+s[2*i+1];
System.out.println(temp[i]);
}
for (int i = 0; i < temp.length; i++) {
for (int i1 = 0; i1 < number; i1++) {
if((temp[i].contains("dn"+(i1+1)+"(t)")&&((i1+1)!=ingex))){
if (temp[i].startsWith("-2")){
temp[i]="-"+temp[i].substring(3,temp[i].length());
} else if (temp[i].startsWith("2")) {
temp[i]=temp[i].substring(2,temp[i].length());
}
}
}
}
for (int i = 0; i < temp.length; i++) {
outa =outa+temp[i];
}
}
return outa;
}
}
注意事项：
比如你想建一个模型，那么请用n1,n2,n3…来表示对应的状态量,dn1,dn2,dn3…来表示对应的速度量。其他物理参数用p1~p9来表示。应该可以理解吧。也就是用这些量描述拉格朗日算子即可。
有不懂的欢迎讨论。
5万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
