matlab的for循环计算多元方程的解,系数变量for循环，求多元一次方程组，运算太慢，怎么优........

本帖最后由 尧遥S 于 2016-2-20 19:47 编辑

我在做一个弹簧力学模型的优化求解，大致是物体质量，质心坐标，弹簧的连接点坐标，弹簧的刚度系数最小值，最大值，步长，物体偏转角等参数由excel文件输入，弹簧的刚度系数为循环变量，求解多元一次方程组，输出满足物体偏转角限制的各组弹簧系数和对应的伸长量。

我的问题简单来说是循环嵌套的矩阵求解(我把解方程组转变为矩阵求解了)速度太慢，怎么优化算法？

我的算法如下：

function chengxu

clear;clc;

s=xlsread('参数.xlsx');%读取参数.xlsx

xa=s(1);ya=s(17);%物质心坐标(xa,ya)

xb=s(2);yb=s(18);%板质心坐标(xb,yb)

x1=s(3);y1=s(19);%弹簧1连接点坐标(x1,y1)

x2=s(4);y2=s(20);%弹簧2连接点坐标(x2,y2)

x3=s(5);y3=s(21);%弹簧3连接点坐标(x3,y3)

x4=s(6);y4=s(22);%弹簧4连接点坐标(x4,y4)

x5=s(7);y5=s(23);%弹簧5连接点坐标(x5,y5)

x6=s(8);y6=s(24);%弹簧6连接点坐标(x6,y6)

x7=s(9);y7=s(25);%弹簧7连接点坐标(x7,y7)

x8=s(10);y8=s(26);%弹簧8连接点坐标(x8,y8)

x9=s(11);y9=s(27);%弹簧9连接点坐标(x9,y9)

Ga=s(49);%物重量

Gb=s(50);%板重量

min=s(13);step=s(14);max=s(15);%弹簧刚度范围及步长min:step:max

theta=s(16);%最大偏转角限制

for k1=min:step:max

for k2=min:step:max

for k3=min:step:max

for k4=min:step:max

for k5=min:step:max

for k6=min:step:max

for k7=min:step:max

for k8=min:step:max

for k9=min:step:max

%物：ΣFz=0,k1*z1+k2*z2+k3*z3+k4*z4+k5*z5=Ga

%物：ΣMx(F)=0,k1*z1*(y1-ya)+(k2*z2+k3*z3)*(y2-ya)+(k4*z4+k5*z5)*(y5-ya)=0

%物：ΣMy(F)=0,(k2*z2+k5*z5)*(x5-xa)+k1*z1*(x1-xa)+(k3*z3+k4*z4)*(x4-xa)=0

%物：∵两条相交直线确定一个平面∴两条对角线各自中点的位置始终重合，即伸长量相等，1/2*(z2+z4)=1/2*(z3+z5)

%物：∵两条相交直线确定一个平面∴两条对角线各自中点的位置始终重合，即伸长量相等，1/2*(1/2*(z3+z4)+z1)=1/2*(z2+z5)

%板：ΣFz=0,k6*z6+k7*z7+k8*z8+k9*z9=k1*z1+k2*z2+k3*z3+k4*z4+k5*z5+Gb,即：k6*z6+k7*z7+k8*z8+k9*z9=Ga+Gb

%板：ΣMx(F)=0,(k6*z6+k7*z7)*(y6-yb)+(k4*z4+k5*z5)*(y5-yb)+(k8*z8+k9*z9)*(y9-yb)+(k2*z2+k3*z3)*(y2-yb)+k1*z1*(y1-yb)=0

%板：ΣMy(F)=0,(k6*z6+k9*z9)*(x6-xb)+(k3*z3+k4*z4)*(x3-xb)+(k7*z7+k8*z8)*(x7-xb)+k1*z1*(x1-xb)+(k2*z2+k5*z5)*(x2-xb)=0

%板：∵两条相交直线确定一个平面∴两条对角线各自中点的位置始终重合，即伸长量相等，1/2*(z6+z8)=1/2*(z7+z9)

%将以上九元一次方程组通过矩阵化简求解(系数矩阵为A，常数项为b，结果为R)

A=[k1 k2 k3 k4 0 0 0 0 0;k1*(y1-ya) k2*(y2-ya) k3*(y2-ya) k4*(y5-ya) k5*(y5-ya) 0 0 0 0;k1*(x1-xa) k2*(x5-xa) k3*(x4-xa) k4*(x4-xa) k5*(x5-xa) 0 0 0 0;0 1 -1 1 -1 0 0 0 0;1 -1 0.5 0.5 -1 0 0 0 0;0 0 0 0 0 k6 k7 k8 k9;k1*(y1-yb) k2*(y2-yb) k3*(y2-yb) k4*(y5-yb) k5*(y5-yb) k6*(y6-yb) k7*(y6-yb) k8*(y9-yb) k9*(y9-yb);k1*(x1-xb) k2*(x2-xb) k3*(x3-xb) k4*(x3-xb) k5*(x2-xb) k6*(x6-xb) k7*(x7-xb) k8*(x7-xb) k9*(x6-xb);0 0 0 0 0 1 -1 1 -1];

b=[Ga;0;0;0;0;Ga+Gb;0;0;0];

R=Ab;

z1=R(1);z2=R(2);z3=R(3);z4=R(4);z5=R(5);z6=R(6);z7=R(7);z8=R(8);z9=R(9);%R=[z1 z2 z3 z4 z5 z6 z7 z8 z9]'

thetax=atand(abs((1/2*(z8+z9)+z5)-(1/2*(z6+z7)+z2))/abs(y5-y2));%此时物x轴方向的偏转角thetax

thetay=atand(abs((1/2*(z7+z8)+1/2*(z3+z4))-(1/2*(z6+z9)+z1))/abs(x3-x1));%此时物y轴方向的偏转角thetay

if (thetax<=theta)&&(thetay

fprintf('k1=%g,k2=%g,k3=%g,k4=%g,k5=%g,k6=%g,k7=%g,k8=%g,k9=%gn',k1,k2,k3,k4,k5,k6,k7,k8,k9)%输出满足条件的各弹簧刚度

fprintf('z1=%g,z2=%g,z3=%g,z4=%g,z5=%g,z6=%g,z7=%g,z8=%g,z9=%gn',z1,z2,z3,z4,z5,z6,z7,z8,z9)%输出此时的各弹簧伸长量

fprintf('thetax=%g,thetay=%gnn',thetax,thetay)%输出此时的x轴方向和y轴方向的偏转角

end

end

end

end

end

end

end

end

end

end

2016-2-20 19:46 上传

点击文件名下载附件

4 KB, 下载次数: 0

2016-2-20 19:46 上传

点击文件名下载附件

9.82 KB, 下载次数: 0