Вывести матрицу пространства состояний из ОДУ в Matlab

Question

Я хотел бы создать коды Matlab, которые автоматически преобразуют ODE в матрицу пространства состояний. Все, что мне нужно, это матрица A и матрица B. Я его частично закончил, и вот как это выглядит. После запуска этого кода я копирую и вставляю результат V, отображаемый в окне команд, в электронную таблицу и вручную нахожу матрицы A и B. Кто-нибудь знает, как я могу создать код, чтобы он автоматически отображал матрицы A и B в командном окне? У меня есть 14 состояний, 4 входа, кстати 4 входа возмущений. Спасибо.

syms x1(t) x2(t) x3(t) x4(t) x5(t) x6(t) x7(t) x8(t) x9(t) x10(t) x11(t) x12(t) x13(t) x14(t) t;
syms Ms Ir Ip m_1 m_2 m_3 m_4 b_1 b_2 b_3 b_4 ks1 ks2 ks3 ks4 a_dis b_dis Tf Tr ku1 ku2 ku3 ku4 zs1 zs2 zs3 zs4;
syms u1(t) u2(t) u3(t) u4(t) w1(t) w2(t) w3(t) w4(t); 

zs1=-Tf*x1-a_dis*x2+x3;
zs2=+Tf*x1-a_dis*x2+x3;
zs3=-Tr*x1+b_dis*x2+x3;
zs4=+Tr*x1+b_dis*x2+x3;

Fd1=b_1*(diff(zs1)-x11);
Fd2=b_2*(diff(zs2)-x12);
Fd3=b_3*(diff(zs3)-x13);
Fd4=b_4*(diff(zs4)-x14);

Fs1=ks1*(zs1-x4);
Fs2=ks2*(zs2-x5);
Fs3=ks3*(zs3-x6);
Fs4=ks4*(zs4-x7);

Fu1= ku1*(x4-w1);
Fu2= ku2*(x5-w2);
Fu3= ku3*(x6-w3);
Fu4= ku4*(x7-w4);

DEq1 = diff(x1)== x8;
DEq2 = diff(x2)== x9;
DEq3 = diff(x3)== x10;
DEq4 = diff(x4)== x11;
DEq5 = diff(x5)== x12;
DEq6 = diff(x6)== x13;
DEq7 = diff(x7)== x14;
DEq8 = Ir*diff(x8)==   -Fd1*Tf+Fd2*Tf-Fd3*Tr+Fd4*Tr-Fs1*Tf+Fs2*Tf-Fs3*Tr+Fs4*Tr+u1*Tf-u2*Tf+u3*Tr-u4*Tr;
DEq9 = Ip*diff(x9)==   -Fd1*a_dis-Fd2*a_dis+Fd3*b_dis+Fd4*b_dis-Fs1*a_dis-Fs2*a_dis+Fs3*b_dis+Fs4*b_dis+u1*a_dis+u2*a_dis-u3*b_dis-u4*b_dis;
DEq10 = Ms*diff(x10)==  -Fd1-Fd2-Fd3-Fd4-Fs1-Fs2-Fs3-Fs4+u1+u2+u3+u4;
DEq11 = m_1*diff(x11)== Fd1+Fs1-Fu1-u1;
DEq12 = m_2*diff(x12)== Fd2+Fs2-Fu2-u2;
DEq13 = m_3*diff(x13)== Fd3+Fs3-Fu3-u3;
DEq14 = m_4*diff(x14)== Fd4+Fs4-Fu4-u4;

[V,S] = odeToVectorField(DEq1,DEq2,DEq3,DEq4,DEq5,DEq6,DEq7,DEq8,DEq9,DEq10,DEq11,DEq12,DEq13,DEq14)

Что отображается в окне команд после выполнения вышеуказанного кода, так это.

V =
 
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       Y[8]
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       Y[9]
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      Y[10]
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      Y[11]
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      Y[12]
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      Y[13]
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      Y[14]
                                                                                     (Tf*u1(t) - Tf*u2(t) + Tr*u3(t) - Tr*u4(t) + Tf^2*b_1*Y[8] + Tf^2*b_2*Y[8] + Tr^2*b_3*Y[8] + Tr^2*b_4*Y[8] + Tf^2*ks1*Y[1] + Tf^2*ks2*Y[1] + Tr^2*ks3*Y[1] + Tr^2*ks4*Y[1] - Tf*b_1*Y[10] + Tf*b_1*Y[11] + Tf*b_2*Y[10] - Tf*b_2*Y[12] - Tr*b_3*Y[10] + Tr*b_4*Y[10] + Tr*b_3*Y[13] - Tr*b_4*Y[14] - Tf*ks1*Y[3] + Tf*ks1*Y[4] + Tf*ks2*Y[3] - Tf*ks2*Y[5] - Tr*ks3*Y[3] + Tr*ks4*Y[3] + Tr*ks3*Y[6] - Tr*ks4*Y[7] + Tf*a_dis*b_1*Y[9] - Tf*a_dis*b_2*Y[9] - Tr*b_3*b_dis*Y[9] + Tr*b_4*b_dis*Y[9] + Tf*a_dis*ks1*Y[2] - Tf*a_dis*ks2*Y[2] - Tr*b_dis*ks3*Y[2] + Tr*b_dis*ks4*Y[2])/Ir
 (a_dis*u1(t) + a_dis*u2(t) - b_dis*u3(t) - b_dis*u4(t) + a_dis^2*b_1*Y[9] + a_dis^2*b_2*Y[9] + b_3*b_dis^2*Y[9] + b_4*b_dis^2*Y[9] + a_dis^2*ks1*Y[2] + a_dis^2*ks2*Y[2] + b_dis^2*ks3*Y[2] + b_dis^2*ks4*Y[2] - a_dis*b_1*Y[10] + a_dis*b_1*Y[11] - a_dis*b_2*Y[10] + a_dis*b_2*Y[12] + b_3*b_dis*Y[10] + b_4*b_dis*Y[10] - b_3*b_dis*Y[13] - b_4*b_dis*Y[14] - a_dis*ks1*Y[3] + a_dis*ks1*Y[4] - a_dis*ks2*Y[3] + a_dis*ks2*Y[5] + b_dis*ks3*Y[3] + b_dis*ks4*Y[3] - b_dis*ks3*Y[6] - b_dis*ks4*Y[7] + Tf*a_dis*b_1*Y[8] - Tf*a_dis*b_2*Y[8] - Tr*b_3*b_dis*Y[8] + Tr*b_4*b_dis*Y[8] + Tf*a_dis*ks1*Y[1] - Tf*a_dis*ks2*Y[1] - Tr*b_dis*ks3*Y[1] + Tr*b_dis*ks4*Y[1])/Ip
                                                                                                                                                                                         (u1(t) + u2(t) + u3(t) + u4(t) - b_1*Y[10] + b_1*Y[11] - b_2*Y[10] - b_3*Y[10] + b_2*Y[12] - b_4*Y[10] + b_3*Y[13] + b_4*Y[14] - ks1*Y[3] + ks1*Y[4] - ks2*Y[3] - ks3*Y[3] + ks2*Y[5] - ks4*Y[3] + ks3*Y[6] + ks4*Y[7] + Tf*b_1*Y[8] - Tf*b_2*Y[8] + Tr*b_3*Y[8] - Tr*b_4*Y[8] + Tf*ks1*Y[1] - Tf*ks2*Y[1] + Tr*ks3*Y[1] - Tr*ks4*Y[1] + a_dis*b_1*Y[9] + a_dis*b_2*Y[9] - b_3*b_dis*Y[9] - b_4*b_dis*Y[9] + a_dis*ks1*Y[2] + a_dis*ks2*Y[2] - b_dis*ks3*Y[2] - b_dis*ks4*Y[2])/Ms
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            -(u1(t) - b_1*Y[10] + b_1*Y[11] - ks1*Y[3] + ks1*Y[4] + ku1*Y[4] - ku1*w1(t) + Tf*b_1*Y[8] + Tf*ks1*Y[1] + a_dis*b_1*Y[9] + a_dis*ks1*Y[2])/m_1
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            -(u2(t) - b_2*Y[10] + b_2*Y[12] - ks2*Y[3] + ks2*Y[5] + ku2*Y[5] - ku2*w2(t) - Tf*b_2*Y[8] - Tf*ks2*Y[1] + a_dis*b_2*Y[9] + a_dis*ks2*Y[2])/m_2
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            -(u3(t) - b_3*Y[10] + b_3*Y[13] - ks3*Y[3] + ks3*Y[6] + ku3*Y[6] - ku3*w3(t) + Tr*b_3*Y[8] + Tr*ks3*Y[1] - b_3*b_dis*Y[9] - b_dis*ks3*Y[2])/m_3
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             (b_4*Y[10] - u4(t) - b_4*Y[14] + ks4*Y[3] - ks4*Y[7] - ku4*Y[7] + ku4*w4(t) + Tr*b_4*Y[8] + Tr*ks4*Y[1] + b_4*b_dis*Y[9] + b_dis*ks4*Y[2])/m_4
 
 

Answer 1

Использование функции Matlab subs - менее уродливый способ, который я нашел для этого. Вы получите матрицу 14x14 sym:

A = repmat(V,1,14);
for j=1:14
    y = double((1:14)==j);
    A(:,j) = subs(A(:,j), 'Y[1]', y(1));
    A(:,j) = subs(A(:,j), 'Y[2]', y(2));
    A(:,j) = subs(A(:,j), 'Y[3]', y(3));
    A(:,j) = subs(A(:,j), 'Y[4]', y(4));
    A(:,j) = subs(A(:,j), 'Y[5]', y(5));
    A(:,j) = subs(A(:,j), 'Y[6]', y(6));
    A(:,j) = subs(A(:,j), 'Y[7]', y(7));
    A(:,j) = subs(A(:,j), 'Y[8]', y(8));
    A(:,j) = subs(A(:,j), 'Y[9]', y(9));
    A(:,j) = subs(A(:,j), 'Y[10]', y(10));
    A(:,j) = subs(A(:,j), 'Y[11]', y(11));
    A(:,j) = subs(A(:,j), 'Y[12]', y(12));
    A(:,j) = subs(A(:,j), 'Y[13]', y(13));
    A(:,j) = subs(A(:,j), 'Y[14]', y(14));
end

Вы можете избавиться от всех строк, используя функцию eval:

A = repmat(V,1,14);
for j=1:14
    y = double((1:14)==j);
    for k=1:14
        A(:,j) = eval(sprintf('subs(A(:,j), ''Y[%d]'', y(%d));',k,k));
    end
end

В столбце j, мы заменяем Y[i] на 1 if i==j else 0.

Test:

>> size(A)

ans =

    14    14

>> A(:,10)
 
ans =
 
                                                                                                           0
                                                                                                           0
                                                                                                           1
                                                                                                           0
                                                                                                           0
                                                                                                           0
                                                                                                           0
                         -(Tf*b_1 - Tf*b_2 + Tr*b_3 - Tr*b_4 - Tf*u1(t) + Tf*u2(t) - Tr*u3(t) + Tr*u4(t))/Ir
 -(a_dis*b_1 + a_dis*b_2 - b_3*b_dis - b_4*b_dis - a_dis*u1(t) - a_dis*u2(t) + b_dis*u3(t) + b_dis*u4(t))/Ip
                                                 -(b_1 + b_2 + b_3 + b_4 - u1(t) - u2(t) - u3(t) - u4(t))/Ms
                                                                               (b_1 - u1(t) + ku1*w1(t))/m_1
                                                                               (b_2 - u2(t) + ku2*w2(t))/m_2
                                                                               (b_3 - u3(t) + ku3*w3(t))/m_3
                                                                               (b_4 - u4(t) + ku4*w4(t))/m_4