前言
建立一个数学模型。通过拉格朗日建模思想可以知道，你只需要确定系统的状态量，并构造拉格朗日算子（实际上，就是算出系统的动能以及势能）即可。
现在，你知道了系统的状态量有五个，也表示出了拉格朗日算子，甚至自己也通过MATLAB计算出了对应的运动方程。但是，你又想通过SIMULINK将模型仿真出来，那没辙，只能根据运动方程挨个找出M、C和G矩阵（向量）的内部元素。
上述工作都做完了之后，你又觉得建模过程太烦了，又要自己不停的替换（至于为什么要替换，去建一遍就知道了），又要自己手动找出矩阵（向量）元素值。你就想了，有没有一个工具，我只需要告诉它系统的状态量个数以及拉格朗日算子表达式，我就可以直接得到最终所有矩阵的元素值呢？
答案是：有的
一、准备工作
该部分为程序的说明，完整程序见文末
1.状态量以及其他参数定义
定义状态量及其速度、加速度项，并定义模型其他参数
syms n1(t) n2(t) n3(t) n4(t) n5(t) %状态量
syms dn1(t) dn2(t) dn3(t) dn4(t) dn5(t) %状态量速度
syms ddn1(t) ddn2(t) ddn3(t) ddn4(t) ddn5(t) %状态量加速度
syms m ml L l D a b g %模型其他参数
2.拉格朗日算子定义
对于该部分的描述，每个模型都不同。但是核心就是根据你前面定义的状态量和系统参数描述出拉格朗日算子。
xp0=0.5*(L*cos(n1)+l*sin(n3)+D-L*cos(n2)-l*sin(n4))-b*cos(n5);%%负载横坐标
yp0=0.5*(L*sin(n1)-l*cos(n3)+L*sin(n2)-l*cos(n4))-b*sin(n5); %%负载纵坐标
dxp0_original=diff(xp0,t); %%对仅包含位置量的表达式求导
dxp0_string=char(dxp0_original); %%将其转化为字符
dxp0=str2sym(replace_diff(dxp0_string)); %%替换预定义的符号变量
dyp0_original=diff(yp0,t); %%对仅包含位置量的表达式求导
dyp0_string=char(dyp0_original); %%将其转化为字符串
dyp0=str2sym(replace_diff(dyp0_string)); %%替换预定义的符号变量
T=(1/6)*ml*L^2*dn1^2+(1/6)*ml*L^2*dn2^2+(1/2)*m*(dxp0^2+dyp0^2)+(1/6)*m*a^2*dn5^2;
P=ml*g*(1/2)*L*sin(n1)+ml*g*(1/2)*L*sin(n2)+m*g*yp0;
L=T-P;%%拉格朗日算子
其中自定义的函数如下：
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(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)');
y1=replace(ans, 'diff(dn5(t), t)', 'ddn5(t)');
end
至此，你得到了系统的拉格朗日算子，下面就让我们开始愉快的建模过程吧！
二、建模
五个状态量，直接甩出五个求解表达式，如下：
eq1=string(simplify(str2sym(replace_diff(char(diff(diff(L,dn1),t))))-diff(L,n1)));
eq2=string(simplify(str2sym(replace_diff(char(diff(diff(L,dn2),t))))-diff(L,n2)));
eq3=string(simplify(str2sym(replace_diff(char(diff(diff(L,dn3),t))))-diff(L,n3)));
eq4=string(simplify(str2sym(replace_diff(char(diff(diff(L,dn4),t))))-diff(L,n4)));
eq5=string(simplify(str2sym(replace_diff(char(diff(diff(L,dn5),t))))-diff(L,n5)));
以eq1为例，string()和simplify()是为了简化公式结构并将其转化为字符串类型，str2sym的作用是将字符串类型的量转化为符号函数，剩下即拉格朗日方程的基本形式。
通过上面的操作，我们得到了完整的系统运动方程。
三、jar包调用
利用JAVA的正则表达式，输入系统运动方程，输出矩阵元素。
编写如下类方法（打包成jar包，供matlab调用）：
package match;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
public class heimaa {
public static void main(String[] args) {
}
public static String[] getVector(String[] equations){
String[][] M = new String[5][5];
String[][] C = new String[5][5];
String[] G = new String[5];
String[] Total = new String[55];
for (int i = 0; i < M.length; i++) {
for (int i1 = 0; i1 < M.length; i1++) {
M[i][i1]="";
}
}
for (int i = 0; i < M.length; i++) {
for (int i1 = 0; i1 < M.length; i1++) {
C[i][i1]="0";
}
}
for (int i = 0; i < G.length; i++) {
G[i]="";
}
for (int j = 0; j < equations.length; j++) {
String str = equations[j].replace(" ","");
str = str.replace("\n","");
Pattern p =Pattern.compile("[+-][(][^\\/\\x22]+[)][/]\\d{1}");
Matcher m = p.matcher(str);
while(m.find()){
String str1 = m.group();
for (int i = 0; i < 5; i++) {
if(str1.contains("*ddn"+(i+1)+"(t)")){
M[j][i]=str1.replace("*ddn"+(i+1)+"(t)","");
str=str.replace(str1,"");
}
if(str1.contains("*dn"+(i+1)+"(t)^2")){
C[j][i]=str1.replace("*dn"+(i+1)+"(t)^2","*dn"+(i+1)+"(t)");
str=str.replace(str1,"");
}
}
}
Pattern p1 =Pattern.compile("[+-][^dd\\x22]+[g][^\\x22]+");
Matcher m1 = p1.matcher(str);
while(m1.find()){
String str3 =m1.group();
G[j]=str3;
str=str.replace(str3,"");
}
for (int i1 = 0; i1 < 5; i1++) {
if(str.contains("*ddn"+(i1+1)+"(t)")){
M[j][i1]=M[j][i1]+"+"+str.replace("*ddn"+(i1+1)+"(t)","");
}
}
}
for (int i = 0; i < 25; i++) {
Total[i]="m"+(i/5+1)+(i%5+1)+"="+M[i/5][i%5];
}
for (int i = 25; i < 50; i++) {
Total[i]="c"+((i-25)/5+1)+((i-25)%5+1)+"="+C[(i-25)/5][(i-25)%5];
}
for (int i = 50; i < 55; i++) {
Total[i]="g"+((i-50)+1)+"="+G[i-50];
}
return Total;
}
}
随后，matlab直接调用jar包，输入上面计算好的运动方程：
obj =match.heimaa;
obj.getVector([eq1,eq2,eq3,eq4,eq5])
下面就愉快的开始仿真吧
四、MATLAB完整程序
%%%%%%首先定义模型状态量及其速度加速度%%%%%%%%
syms n1(t) n2(t) n3(t) n4(t) n5(t) %状态量
syms dn1(t) dn2(t) dn3(t) dn4(t) dn5(t) %状态量速度
syms ddn1(t) ddn2(t) ddn3(t) ddn4(t) ddn5(t) %状态量加速度
syms m ml L l D a b g %m:DMB质量 ml:悬臂质量 L:悬臂长度 l:悬绳长度 D:悬臂底端点水平距离 a:DMB长度 b:0.5*DMB宽度
%%%%%%负载的坐标与速度%%%%%%
xp0=0.5*(L*cos(n1)+l*sin(n3)+D-L*cos(n2)-l*sin(n4))-b*cos(n5);%%负载横坐标
yp0=0.5*(L*sin(n1)-l*cos(n3)+L*sin(n2)-l*cos(n4))-b*sin(n5); %%负载纵坐标
dxp0_original=diff(xp0,t);%%对仅包含位置量的表达式求导
dxp0_string=char(dxp0_original);%%将其转化为字符
dxp0=str2sym(replace_diff(dxp0_string));%%替换预定义的符号变量
dyp0_original=diff(yp0,t);%%对仅包含位置量的表达式求导
dyp0_string=char(dyp0_original);%%将其转化为字符串
dyp0=str2sym(replace_diff(dyp0_string));%%替换预定义的符号变量
%%%%%%动能与使能(直接写出来,不涉及求导)%%%%%%
T=(1/6)*ml*L^2*dn1^2+(1/6)*ml*L^2*dn2^2+(1/2)*m*(dxp0^2+dyp0^2)+(1/6)*m*a^2*dn5^2;
P=ml*g*(1/2)*L*sin(n1)+ml*g*(1/2)*L*sin(n2)+m*g*yp0;
L=T-P;%%拉格朗日算子
%%%%%方程计算%%%%%
eq1=string(simplify(str2sym(replace_diff(char(diff(diff(L,dn1),t))))-diff(L,n1)));
eq2=string(simplify(str2sym(replace_diff(char(diff(diff(L,dn2),t))))-diff(L,n2)));
eq3=string(simplify(str2sym(replace_diff(char(diff(diff(L,dn3),t))))-diff(L,n3)));
eq4=string(simplify(str2sym(replace_diff(char(diff(diff(L,dn4),t))))-diff(L,n4)));
eq5=string(simplify(str2sym(replace_diff(char(diff(diff(L,dn5),t))))-diff(L,n5)));
%%生成矩阵元素%%%
obj =match.heimaa;
obj.getVector([eq1,eq2,eq3,eq4,eq5])
运行结果：
'm11=+(L^2*ml)/3+(L^2*m)/4'
'm12=+(L^2*m*cos(n1(t)+n2(t)))/4'
'm13=-(L*l*m*sin(n1(t)-n3(t)))/4'
'm14=+(L*l*m*sin(n1(t)+n4(t)))/4'
'm15=-(L*b*m*cos(n1(t)-n5(t)))/2'
'm21=+(L^2*m*cos(n1(t)+n2(t)))/4'
'm22=+(L^2*ml)/3+(L^2*m)/4'
'm23=+(L*l*m*sin(n2(t)+n3(t)))/4'
'm24=-(L*l*m*sin(n2(t)-n4(t)))/4'
'm25=-(L*b*m*cos(n2(t)+n5(t)))/2'
'm31=-(L*l*m*sin(n1(t)-n3(t)))/4'
'm32=+(L*l*m*sin(n2(t)+n3(t)))/4'
'm33=+(l^2*m)/4'
'm34=-(l^2*m*cos(n3(t)+n4(t)))/4'
'm35=-(b*l*m*sin(n3(t)-n5(t)))/2'
'm41=+(L*l*m*sin(n1(t)+n4(t)))/4'
'm42=-(L*l*m*sin(n2(t)-n4(t)))/4'
'm43=-(l^2*m*cos(n3(t)+n4(t)))/4'
'm44=+(l^2*m)/4'
'm45=-(b*l*m*sin(n4(t)+n5(t)))/2'
'm51=-(L*b*m*cos(n1(t)-n5(t)))/2'
'm52=-(L*b*m*cos(n2(t)+n5(t)))/2'
'm53=-(b*l*m*sin(n3(t)-n5(t)))/2'
'm54=-(b*l*m*sin(n4(t)+n5(t)))/2'
'm55=+(a^2*m)/3+b^2*m'
'c11=0'
'c12=-(L^2*m*sin(n1(t)+n2(t))*dn2(t))/4'
'c13=+(L*l*m*dn3(t)*cos(n1(t)-n3(t)))/4'
'c14=+(L*l*m*cos(n1(t)+n4(t))*dn4(t))/4'
'c15=-(L*b*m*dn5(t)*sin(n1(t)-n5(t)))/2'
'c21=-(L^2*m*sin(n1(t)+n2(t))*dn1(t))/4'
'c22=0'
'c23=+(L*l*m*cos(n2(t)+n3(t))*dn3(t))/4'
'c24=+(L*l*m*dn4(t)*cos(n2(t)-n4(t)))/4'
'c25=+(L*b*m*sin(n2(t)+n5(t))*dn5(t))/2'
'c31=-(L*l*m*dn1(t)*cos(n1(t)-n3(t)))/4'
'c32=+(L*l*m*cos(n2(t)+n3(t))*dn2(t))/4'
'c33=0'
'c34=+(l^2*m*sin(n3(t)+n4(t))*dn4(t))/4'
'c35=+(b*l*m*dn5(t)*cos(n3(t)-n5(t)))/2'
'c41=+(L*l*m*cos(n1(t)+n4(t))*dn1(t))/4'
'c42=-(L*l*m*dn2(t)*cos(n2(t)-n4(t)))/4'
'c43=+(l^2*m*sin(n3(t)+n4(t))*dn3(t))/4'
'c44=0'
'c45=-(b*l*m*cos(n4(t)+n5(t))*dn5(t))/2'
'c51=+(L*b*m*dn1(t)*sin(n1(t)-n5(t)))/2'
'c52=+(L*b*m*sin(n2(t)+n5(t))*dn2(t))/2'
'c53=-(b*l*m*dn3(t)*cos(n3(t)-n5(t)))/2'
'c54=-(b*l*m*cos(n4(t)+n5(t))*dn4(t))/2'
'c55=0'
'g1=+(L*g*m*cos(n1(t)))/2+(L*g*ml*cos(n1(t)))/2'
'g2=+(L*g*m*cos(n2(t)))/2+(L*g*ml*cos(n2(t)))/2'
'g3=+(g*l*m*sin(n3(t)))/2'
'g4=+(g*l*m*sin(n4(t)))/2'
'g5=-b*g*m*cos(n5(t))'
总结
本文涉及到matlab调用jar包、JAVA正则表达式、拉格朗日建模等知识点，值得注意的是，JAVA编写的程序可能不适用于某些模型，但总体思想一致，欢迎各位对新人的第一篇文章进行指导。
1065
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
