1. Tuyển Mod quản lý diễn đàn. Các thành viên xem chi tiết tại đây

SOS! PTVP trong matlab. Bất kì ai trả lời cũng được vote

Chủ đề trong 'Cơ khí - Tự động hoá' bởi chutmayman, 25/10/2005.

  1. 0 người đang xem box này (Thành viên: 0, Khách: 0)
  1. chutmayman

    chutmayman Thành viên mới

    Tham gia ngày:
    12/10/2005
    Bài viết:
    1.033
    Đã được thích:
    1
    SOS! PTVP trong matlab. Bất kì ai trả lời cũng được vote

    E đang có chút thắc mắc về giải hệ ptvp phi tuyến trong matlab.
    Các bác giúp e nhé.
    Bác nào giúp e, e sẽ gửi mail mfile đến.
    Xin cảm ơn và hậu tạ vote5*
  2. WJT

    WJT Thành viên mới

    Tham gia ngày:
    09/10/2005
    Bài viết:
    492
    Đã được thích:
    4
    Bạn có thắc mắc gì thì cứ post lên đây mới có nhiều chuyên gia để ý chứ!
    WJT.
  3. chutmayman

    chutmayman Thành viên mới

    Tham gia ngày:
    12/10/2005
    Bài viết:
    1.033
    Đã được thích:
    1
    E gửi các bác 2 file mfile. Các bác xem giúp nhé.
    Khi e chạy file pl1 thi nó báo lỗi như sau:
    ??? Input argument ''t'' is undefined.
    Error in ==> D:MATLAB6p5workQTQD.m
    On line 7 ==> uA=sqrt(2)*Um*sin(2*pi*fv*t);
    Error in ==> D:MATLAB6p5workPL1.M
    On line 37 ==> [t,si] = ode45(QTQD,[0 0.7],[0;0;0;0;0]);
  4. chutmayman

    chutmayman Thành viên mới

    Tham gia ngày:
    12/10/2005
    Bài viết:
    1.033
    Đã được thích:
    1
    Đây là file chính
    %clc
    %clear
    global r1 F1 x10 x20 Xo ktg Dx2s Omeo r20 r21 is p G ...
    a1 a2 b1 b2 s X1 X2 r2 x2s Mc Um fv;
    a1 = 0; a2 = 0; b1 = 0; b2 = 0;
    is = 0; f = 50;
    x10 = 0.98; %so lieu dong co
    x20 = 0.92; %so lieu dong co
    x2s = 0.81; %so lieu dong co
    x2bh = 0.755; %so lieu dong co
    r20 = 0.235; %so lieu dong co
    r21 = 0.376; %so lieu dong co
    r1 = 0.4; %so lieu dong co
    Xo = 22; %so lieu dong co
    p = 2; %so lieu dong co,so doi cuc
    J = 0.125; %so lieu dong co, momen quan tinh
    ktg = 1/115;
    Dx2s = x20 - x2s;
    Omeo = 314; %Toc do goc dong bo omega_0
    F1 = p/J;
    fv = 50; %Tan so can thay doi
    %----------------------------------------------------------------------------------
    %U/f khong doi, Mc = const
    ku = 220/f;
    Um = ku*fv;
    %----------------------------------------------------------------------------------
    %U/f^2 khong doi, Mc ~ Omega^2
    %ku = 220/f^2;
    %Um = ku*fv^2;
    %----------------------------------------------------------------------------------
    %U/sqrt(f^3) khong doi, Mc ~ Omega
    %ku = 220/sqrt(f^3);
    %Um = ku*sqrt(fv^3);
    %----------------------------------------------------------------------------------
    [t,si] = ode45(QTQD,[0 0.7],[0;0;0;0;0]);
    uA=sqrt(2)*Um*sin(2*pi*fv*t'');
    uB=sqrt(2)*Um*sin(2*pi*fv*t''+2*pi/3);
    uC=sqrt(2)*Um*sin(2*pi*fv*t''-2*pi/3);
    u1_anpha = uA;
    u1_beta = (sqrt(3)/3)*(uC-uB);
    i1_anpha = si:-),1)*a1 - si:-),3)*a2;
    i2_anpha = si:-),3)*b1 - si:-),1)*b2;
    i1_beta = si:-),2)*a1 - si:-),4)*a2;
    i2_beta = si:-),4)*b1 - si:-),2)*b2;
    M = G*(si:-),3).*i1_beta - si:-),4).*i1_anpha);
    Di1_anpha = a1*u1_anpha - a1*r1*i1_anpha + a2*r2*i2_anpha + a2*si:-),5).*si:-),4);
    Di1_beta = a1*u1_beta - a1*r1*i1_beta + a2*r2*i2_beta - a2*si:-),5).*si:-),3);
    Di2_anpha = -b2*(u1_anpha - r1*i1_anpha) - b1*(r2*i2_anpha + si:-),5).*si:-),4));
    Di2_beta = -b1*(u1_beta-r1*i1_beta) + b1*(-r2*i2_beta + si:-),5).*si:-),3));
    pd1 = r1*(i1_anpha.^2 + i1_beta.^2);
    pd2 = r2*(i2_anpha.^2 + i2_beta.^2);
    pd = pd1 + pd2;
    im_anpha = i1_anpha + i2_anpha;
    im_beta = i1_beta + i2_beta;
    Dim_anpha = Di1_anpha + Di2_anpha;
    Dim_beta = Di1_beta + Di2_beta;
    pdg1 = (X1/Omeo)*(i1_anpha.*Di1_anpha + i1_beta.*Di1_beta);
    pdg2 = (X2/Omeo)*(i2_anpha.*Di2_anpha + i2_beta.*Di2_beta);
    pdg = pdg1 + pdg2;
    pm = (Xo/Omeo)*(im_anpha.*Dim_anpha + im_beta.*Dim_beta);
    P1 = u1_anpha.*i1_anpha+u1_beta.*i1_beta;
    Tp = pd + pdg + pm;
    temp4 = size(t);
    kt = temp4(1,1);
    %Tinh thoi gian tich phan
    Count = 0;
    for i=1:kt
    if si(i,5) < 0.95*si(kt,5);
    Count = Count + 1;
    else
    break
    end
    end

    %Tinh dien nang tieu thu
    for i=1:kt
    if t(i) < t(count)
    temp(i) = t(i);
    end
    end
    temp2 = size(temp);
    Dtp = temp2(1,2);
    A1 = trapz(t(1:Dtp),P1(1:Dtp));
    Adg = trapz(t(1:Dtp),pdg(1:Dtp));
    ATp = trapz(t(1:Dtp),Tp(1:Dtp));
    Ad = trapz(t(1:Dtp),pd(1:Dtp));
    Phan_tram_Adg_A1 = Adg/A1;
    Phan_tram_Ad_A1 = Ad/A1;
    Phan_tram_Adg_ATp = Adg/ATp;
    %----------------------------------------------------------------------------------
    figure(1);
    plot(t,si:-),5),t,M);
    hold on;
    title([''Thoi diem dong co o che do xac lap: '' num2str(t(Dtp))]);
    xlabel(''thoi gian'');
    ylabel(''MOMEN va TOC DO'');
    %legend(''Tocdo'',''Momen'');
    grid on;
    figure(2);
    subplot(2,1,1);
    plot(t,u1_anpha);
    hold on;
    title(''Dien ap pha thuc'');
    xlabel(''Thoi gian (giay)'');
    ylabel(''Dien ap (V)'');
    grid on;
    subplot(2,1,2);
    plot(t,i1_anpha);
    hold on;
    title(''Dong dien pha A stato'');
    xlabel(''Thoi gian (giay)'');
    ylabel(''Dong dien (A)'');
    grid on;
    figure(3);
    subplot(3,1,1);
    plot(t, pdg);
    hold on;
    title(''Ton hao qua do'');
    xlabel(''Thoi gian (giay)'');
    ylabel(''Cong suat (W)'');
    grid on;
    subplot(3,1,2);
    plot(t, pd);
    hold on;
    title(''Ton hao dong trong qua trinh qua do'');
    xlabel(''Thoi gian (giay)'');
    ylabel(''Cong suat (W)'');
    grid on;
    subplot(3,1,3);
    plot(t, Tp);
    hold on;
    title(''Tong ton hao trong qua trinh qua do'');
    xlabel(''Thoi gian (giay)'');
    ylabel(''Cong suat (W)'');
    grid on;
    Đây là chương trình con
    unction dsi = QTQD(t,si);
    global r1 F1 i1_anpha i1_beta x10 Xo ktg x20 Dx2s Omeo r20 r21 is p G ...
    a1 a2 b1 b2 X1 X2 r2 s x2s X1 X2 Mc Um fv;

    % Tinh cac tham so
    uA=sqrt(2)*Um*sin(2*pi*fv*t);
    uB=sqrt(2)*Um*sin(2*pi*fv*t+2*pi/3);
    uC=sqrt(2)*Um*sin(2*pi*fv*t-2*pi/3);
    u1_anpha = uA;
    u1_beta = (sqrt(3)/3)*(uC-uB);
    % Tinh toan phan tham so dong
    s = (1/Omeo)*(Omeo - si(5));
    is = sqrt((si(1)*a1 - si(3)*a2)^2 + (si(2)*a1 - si(4)*a2)^2);
    X1 = x10 - 1.15*(1e-2)*ktg*is - 0.25;
    X2 = x20 - 0.11*(1e-2)*ktg*is - s*Dx2s;
    r2 = r20 + s*(r21 - r20);
    %--------------------------------------------------------------------------
    Xs = Xo + X1;
    Xr = Xo + X2;
    a1 = (Omeo*Xr) / (Xs*Xr - Xo^2);
    a2 = a1*(Xo/Xr);
    b1 = (Omeo*Xs) / (Xs*Xr - Xo^2);
    b2 = b1*(Xo/Xs);
    kr = Xo/Xr;
    G = (3*p*kr)/2;
    % Xac dinh phu tai can tinh
    Mc = 0;
    %Mc = 92.2;
    %Mc = 0.00093*si(5)^2;
    %Mc = 0.2936*si(5);
    %--------------------------------------------------------------------------
    i1_anpha = si(1)*a1 - si(3)*a2;
    i2_anpha = si(3)*b1 - si(1)*b2;
    i1_beta = si(2)*a1 - si(4)*a2;
    i2_beta = si(4)*b1 - si(2)*b2;
    M = G*(si(3).*(si(2)*a1 - si(4)*a2) - si(4).*(si(1)*a1 - si(3)*a2));
    dsi = zeros(5,1);
    % Xac dinh he ham vi phan
    dsi(1) = u1_anpha - r1*i1_anpha;
    dsi(2) = u1_beta - r1*i1_beta;
    dsi(3) = -r2*i2_anpha - si(5)*si(4);
    dsi(4) = -r2*i2_beta - si(5)*si(3);
    dsi(5) = F1*(M-Mc);
    Được chutmayman sửa chữa / chuyển vào 21:23 ngày 26/10/2005
  5. chutmayman

    chutmayman Thành viên mới

    Tham gia ngày:
    12/10/2005
    Bài viết:
    1.033
    Đã được thích:
    1
    Các chuyên gia giúp em với.
  6. LE_ROB

    LE_ROB Thành viên mới

    Tham gia ngày:
    22/08/2004
    Bài viết:
    220
    Đã được thích:
    0
    CT của bạn bị lặp va số loi dùng biến thế nên không chạy được là phải.
    VD:
    - Lời gọi của hàm ode45 tham chiếu tới hàm QTQD, trong hàm này lại dùng biến u1_alpha, biến được tính phía sau lời gọi hàm ode45.
    - Trong lời gọi của QTQD cần có tham số cho nó (cái này do bạn đã định nghĩa cho QTQD).
    Tớ copy lại và trình bầy qua CT của bạn. Tuy nhiên nó vẫn chưa chạy được đâu, bạn xem lại nhé.
    function[] = un()
    global r1 F1 x10 x20 Xo ktg Dx2s Omeo r20 r21 is p G ...
    a1 a2 b1 b2 s X1 X2 r2 x2s Mc Um fv;
    a1 = 0; a2 = 0; b1 = 0; b2 = 0;
    is = 0; f = 50;
    x10 = 0.98; %so lieu dong co
    x20 = 0.92; %so lieu dong co
    x2s = 0.81; %so lieu dong co
    x2bh = 0.755; %so lieu dong co
    r20 = 0.235; %so lieu dong co
    r21 = 0.376; %so lieu dong co
    r1 = 0.4; %so lieu dong co
    Xo = 22; %so lieu dong co
    p = 2; %so lieu dong co,so doi cuc
    J = 0.125; %so lieu dong co, momen quan tinh
    ktg = 1/115;
    Dx2s = x20 - x2s;
    Omeo = 314; %Toc do goc dong bo omega_0
    F1 = p/J;
    fv = 50; %Tan so can thay doi
    %----------------------------------------------------------------------------------
    %U/f khong doi, Mc = const
    ku = 220/f;
    Um = ku*fv;
    %----------------------------------------------------------------------------------
    %U/f^2 khong doi, Mc ~ Omega^2
    %ku = 220/f^2;
    %Um = ku*fv^2;
    %----------------------------------------------------------------------------------
    %U/sqrt(f^3) khong doi, Mc ~ Omega
    %ku = 220/sqrt(f^3);
    %Um = ku*sqrt(fv^3);
    %----------------------------------------------------------------------------------
    [t,si] = ode45(QTQD,[0, 0.7],[0,0,0,0,0]); %%%%% THAM SO CUA QTQD
    uA=sqrt(2)*Um*sin(2*pi*fv*t'''');
    uB=sqrt(2)*Um*sin(2*pi*fv*t''''+2*pi/3);
    uC=sqrt(2)*Um*sin(2*pi*fv*t''''-2*pi/3);
    u1_anpha = uA;
    u1_beta = (sqrt(3)/3)*(uC-uB);
    i1_anpha = si:-),1)*a1 - si:-),3)*a2;
    i2_anpha = si:-),3)*b1 - si:-),1)*b2;
    i1_beta = si:-),2)*a1 - si:-),4)*a2;
    i2_beta = si:-),4)*b1 - si:-),2)*b2;
    M = G*(si:-),3).*i1_beta - si:-),4).*i1_anpha);
    Di1_anpha = a1*u1_anpha - a1*r1*i1_anpha + a2*r2*i2_anpha + a2*si:-),5).*si:-),4);
    Di1_beta = a1*u1_beta - a1*r1*i1_beta + a2*r2*i2_beta - a2*si:-),5).*si:-),3);
    Di2_anpha = -b2*(u1_anpha - r1*i1_anpha) - b1*(r2*i2_anpha + si:-),5).*si:-),4));
    Di2_beta = -b1*(u1_beta-r1*i1_beta) + b1*(-r2*i2_beta + si:-),5).*si:-),3));
    pd1 = r1*(i1_anpha.^2 + i1_beta.^2);
    pd2 = r2*(i2_anpha.^2 + i2_beta.^2);
    pd = pd1 + pd2;
    im_anpha = i1_anpha + i2_anpha;
    im_beta = i1_beta + i2_beta;
    Dim_anpha = Di1_anpha + Di2_anpha;
    Dim_beta = Di1_beta + Di2_beta;
    pdg1 = (X1/Omeo)*(i1_anpha.*Di1_anpha + i1_beta.*Di1_beta);
    pdg2 = (X2/Omeo)*(i2_anpha.*Di2_anpha + i2_beta.*Di2_beta);
    pdg = pdg1 + pdg2;
    pm = (Xo/Omeo)*(im_anpha.*Dim_anpha + im_beta.*Dim_beta);
    P1 = u1_anpha.*i1_anpha+u1_beta.*i1_beta;
    Tp = pd + pdg + pm;
    temp4 = size(t);
    kt = temp4(1,1);
    %Tinh thoi gian tich phan
    Count = 0;
    for i=1:kt
    if si(i,5) < 0.95*si(kt,5);
    Count = Count + 1;
    else
    break
    end
    end
    %Tinh dien nang tieu thu
    for i=1:kt
    if t(i) < t(count)
    temp(i) = t(i);
    end
    end
    temp2 = size(temp);
    Dtp = temp2(1,2);
    A1 = trapz(t(1:Dtp),P1(1:Dtp));
    Adg = trapz(t(1:Dtp),pdg(1:Dtp));
    ATp = trapz(t(1:Dtp),Tp(1:Dtp));
    Ad = trapz(t(1:Dtp),pd(1:Dtp));
    Phan_tram_Adg_A1 = Adg/A1;
    Phan_tram_Ad_A1 = Ad/A1;
    Phan_tram_Adg_ATp = Adg/ATp;
    %----------------------------------------------------------------------------------
    figure(1);
    plot(t,si:-),5),t,M);
    hold on;
    title([''''Thoi diem dong co o che do xac lap: '''' num2str(t(Dtp))]);
    xlabel(''''thoi gian'''');
    ylabel(''''MOMEN va TOC DO'''');
    %legend(''''Tocdo'''',''''Momen'''');
    grid on;
    figure(2);
    subplot(2,1,1);
    plot(t,u1_anpha);
    hold on;
    title(''''Dien ap pha thuc'''');
    xlabel(''''Thoi gian (giay)'''');
    ylabel(''''Dien ap (V)'''');
    grid on;
    subplot(2,1,2);
    plot(t,i1_anpha);
    hold on;
    title(''''Dong dien pha A stato'''');
    xlabel(''''Thoi gian (giay)'''');
    ylabel(''''Dong dien (A)'''');
    grid on;
    figure(3);
    subplot(3,1,1);
    plot(t, pdg);
    hold on;
    title(''''Ton hao qua do'''');
    xlabel(''''Thoi gian (giay)'''');
    ylabel(''''Cong suat (W)'''');
    grid on;
    subplot(3,1,2);
    plot(t, pd);
    hold on;
    title(''''Ton hao dong trong qua trinh qua do'''');
    xlabel(''''Thoi gian (giay)'''');
    ylabel(''''Cong suat (W)'''');
    grid on;
    subplot(3,1,3);
    plot(t, Tp);
    hold on;
    title(''''Tong ton hao trong qua trinh qua do'''');
    xlabel(''''Thoi gian (giay)'''');
    ylabel(''''Cong suat (W)'''');
    grid on;
    %??y l? chuong tr?nh con
    function[dsi] = QTQD(t,si);
    global r1 F1 i1_anpha i1_beta x10 Xo ktg x20 Dx2s Omeo r20 r21 is p G a1 a2 b1 b2 X1 X2 r2 s x2s X1 X2 Mc Um fv;
    % Tinh cac tham so
    uA=sqrt(2)*Um*sin(2*pi*fv*t);
    uB=sqrt(2)*Um*sin(2*pi*fv*t+2*pi/3);
    uC=sqrt(2)*Um*sin(2*pi*fv*t-2*pi/3);
    u1_anpha = uA;
    u1_beta = (sqrt(3)/3)*(uC-uB);
    % Tinh toan phan tham so dong
    s = (1/Omeo)*(Omeo - si(5));
    is = sqrt((si(1)*a1 - si(3)*a2)^2 + (si(2)*a1 - si(4)*a2)^2);
    X1 = x10 - 1.15*(1e-2)*ktg*is - 0.25;
    X2 = x20 - 0.11*(1e-2)*ktg*is - s*Dx2s;
    r2 = r20 + s*(r21 - r20);
    %--------------------------------------------------------------------------
    Xs = Xo + X1;
    Xr = Xo + X2;
    a1 = (Omeo*Xr) / (Xs*Xr - Xo^2);
    a2 = a1*(Xo/Xr);
    b1 = (Omeo*Xs) / (Xs*Xr - Xo^2);
    b2 = b1*(Xo/Xs);
    kr = Xo/Xr;
    G = (3*p*kr)/2;
    % Xac dinh phu tai can tinh
    Mc = 0;
    %Mc = 92.2;
    %Mc = 0.00093*si(5)^2;
    %Mc = 0.2936*si(5);
    %--------------------------------------------------------------------------
    i1_anpha = si(1)*a1 - si(3)*a2;
    i2_anpha = si(3)*b1 - si(1)*b2;
    i1_beta = si(2)*a1 - si(4)*a2;
    i2_beta = si(4)*b1 - si(2)*b2;
    M = G*(si(3).*(si(2)*a1 - si(4)*a2) - si(4).*(si(1)*a1 - si(3)*a2));
    dsi = zeros(5,1);
    % Xac dinh he ham vi phan
    dsi(1) = u1_anpha - r1*i1_anpha;
    dsi(2) = u1_beta - r1*i1_beta;
    dsi(3) = -r2*i2_anpha - si(5)*si(4);
    dsi(4) = -r2*i2_beta - si(5)*si(3);
    dsi(5) = F1*(M-Mc);
    Được le_rob sửa chữa / chuyển vào 01:41 ngày 02/11/2005
  7. chutmayman

    chutmayman Thành viên mới

    Tham gia ngày:
    12/10/2005
    Bài viết:
    1.033
    Đã được thích:
    1
    Cảm ơn nhiều nhé.
    Bạn có thể cho địa chỉ liên hệ để mình trao đổi trực tiếp k?
    Mà bạn ở đâu nhỉ?
  8. LE_ROB

    LE_ROB Thành viên mới

    Tham gia ngày:
    22/08/2004
    Bài viết:
    220
    Đã được thích:
    0
    Hi, tớ không ở VN có gì bạn cứ PS trong TTVNOL cho tớ, thỉnh thoảng tớ vào đây chơi mà.
  9. LE_ROB

    LE_ROB Thành viên mới

    Tham gia ngày:
    22/08/2004
    Bài viết:
    220
    Đã được thích:
    0
    [/quote]
    Hi, tớ không ở VN có gì bạn cứ PS trong TTVNOL cho tớ, thỉnh thoảng tớ vào đây chơi mà.
    [/quote]
    Hoặc bạn gửi cho tớ tới hộp thư: let_le_rob@yahoo.com

Chia sẻ trang này