›
›

Bài giảng Tự động hóa tính toán thiết kế tàu

pdf 173 trang phuongnguyen 2900
Download
Bạn đang xem 20 trang mẫu của tài liệu "Bài giảng Tự động hóa tính toán thiết kế tàu", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

  • pdfbai_giang_tu_dong_hoa_tinh_toan_thiet_ke_tau.pdf

Nội dung text: Bài giảng Tự động hóa tính toán thiết kế tàu

  1. TRẦN CÔNG NGHỊ TỰ ĐỘNG HÓA TÍNH TOÁN THIẾT KẾ TÀU ĐẠI HỌC GIAO THÔNG VẬN TẢI TP HỒ CHÍ MINH 7 - 2009
  2. Trang để trống
  3. Trần công nghị TỰ ĐỘNG HÓA TÍNH TOÁN THIẾT KẾ TÀU THÀNH PHỐ HỒ CHÍ MINH 2009 ĐẠI HỌC GIAO THÔNG VẬN TẢI TP HỒ CHÍ MINH
  4. Mục lục Mở đầu 5 Chương 1: Phương pháp tính và tự động hóa ính toán thiết kế tàu 6 1.1 Nội suy Lagrange 6 1.2 Tích phân một lớp 8 1.3 Đa thức Legendre 11 1.4 Đa thức Tchebyshev 14 1.5 Tìm nghiệm bằng phương pháp chia đôi đoạn có nghiệm 15 1.6 Phương pháp tổng nhỏ nhất các bình phương 16 1.7 Qui hoạch tuyến tính 21 1.8 Qui hoạch phi tuyến 27 1.8.1 Hàm một biến 28 1.8.2 Hàm nhiều biến 29 1.8.3 Xác định min/max hàm một biến 30 1.8.4 Phương pháp sử dụng gradient 34 1.8.5 Phương pháp tìm trực tiếp (không qua giai đoạn tính gradient). 39 1.8.6 Phương pháp dùng hàm phạt penalty 43 Chương 2: Tính nổi và tính ổn định tàu 51 2.1 Tính nổi tàu thủy 51 2.1.1 Kích thước chính và các hệ số thân tàu 51 2.1.2 Tỷ lệ Bonjean 55 2.1.3 Tính thể tích phân chìm và cá đại lượng liên quan thể tích 55 2.1.4 Tính các đường thủy tĩnh trên máy cá nhân 58 2.1.5 Biểu đồ Firsov 60 2.2 Ổn định tàu. 61 2.2.1 Ổn định ngang ban đầu. 61 2.2.2 Ổn định tại góc nghiêng lớn. 62 2.2.3 Đồ thị ổn định. 65 2.3 Thuật toán xác lập họ đường Cross Curves (pantokaren) 65 2.4 Giới thiệu chương trình tính tính nổi tàu thủy 70 Chương 3: Sức cản vỏ tàu 78 3.1 Sức cản vỏ tàu 77 3.2 Công suất hữu hiệu 81 3.3 Các phương pháp kinh nghiệm tính sức cản vỏ tàu 81 Chương 4: Thiết kế chân vịt tàu thủy 95 4.1 Đặc tính hình học chân vịt 95 4.2 Vẽ chân vịt 97 4.3 Đặc tính thủy động lực 98 4.4 Đồ thị thiết kế chân vịt 102 4.5 Tính hệ số dòng theo, hệ số lực hút 106 4.6 Xâm thực chân vịt 109 4.7 Độ bền cánh chân vịt 114 4.8 Thiết kế chân vịt bước cố định 118 4.9 Lập chương trình thiết kế chân vịt tàu 126 4.10 Vẽ chân vịt trên máy PC 135 3
  5. Chương 5: Thiết kế tối ưu tàu thủy 148 5.1 Đánh giá các chỉ tiêu kinh tế – kỹ thuật của tàu 148 5.2 Sơ đồ tính hiệu quả kinh tế 150 5.3 Tự động thiết kế tàu vận tải 151 Tài liệu tham khảo 172 4
  6. Mở đầu “Tự động hóa tính toán, thiết kế tàu” trình bày cách tính toán, thuật toán phục vụ việc lập chương trình tính tính năng tàu thủy, tính di chuyển, thiết bị đẩy tàu và tự động hóa vẽ tàu. Sau mười năm sử dụng sách cho chuyên đề này những người viết chỉnh, sửa, viết lại phù hợp thực tế. Sửa chữa và bổ sung lần này nhằm làm cho tài liệu sát đề cương giảng dạy và học tập tại trường Đại học Giao thông Vận tải Tp Hồ Chí Minh. Hy vọng rằng sách có ích cho những người đang theo học đóng tàu và công trình nổi cũng như các đồng nghiệp đang làm việc trong cùng lĩnh vực. Thành phố Hồ Chí Minh tháng 6 năm 2009. Người viết “Mở đầu” lần in thứ nhất “Tự động hóa tính toán, thiết kế và đóng tàu” bao gồm hướng dẫn tính toán, chương trình tính phục vụ những môn học tàu thủy tại trường đại học. Những đề tài trong tài liệu này: Thiết kế tàu, Tính nổi và tính ổn định, Sức cản vỏ tàu và thiết bị đẩy tàu, Qui hoạch tuyến tính, qui hoạch phi tuyến và ứng dụng của lý thuyết này vào thiết kế tối ưu tàu thủy, Spline và ứng dụng trong vẽ đường hình, khai triển vỏ tàu. Tài liệu được bố trí theo cách tiện lợi cho người đọc. Mở đầu mỗi chương bạn đọc có điều kiện ôn lại những hiểu biết cần thiết về các phương pháp tính liên quan đến nội dung của chương, có điều kiện làm quen chương trình tính viết bằng ngôn ngữ C áp dụng trong tính toán. Các chương trình nhỏ này còn được dùng cho những vấn đề liên quan với ngành tàu. Nội dung mỗi chương chỉ gồm những kiến thức đã được truyền đạt trong trường đại học chuyên ngành. Trên cơ sở những vấn đề đang được trình bày bạn đọc tìm hiểu thêm giải thuật xử lý những bài toán cụ thể đang đặt ra và cách hoàn thiện một chương trình máy tính dựa vào giải thuật vừa có. Tài liệu có thể giúp ích cho sinh viên khoa đóng tàu, kỹ sư làm việc trong lĩnh vực đóng sửa tàu, thiết kế, nghiên cứu tàu cùng đông đảo bạn đọc quan tâm đến tàu thủy khi tính toán tính năng, như tính nổi, ổn định, tính sức cản, chọn máy phù hợp, thiết kế mới, lập phương án đóng mới, lập phương án sửa chữa tàu vv Trong quá trình biên soạn tài liệu những người làm công tác chuẩn bị nhận được sự giúp đỡ chân tình và thiết thực của ban giám hiệu phân hiệu Đại học Hàng hải, phân khoa đóng tàu, bạn cùng nghề và những bạn bè xa, gần. Những đóng góp quí giá về nội dung, về biên soạn và hiệu chỉnh tài liệu, hiệu chỉnh các bản in thử vv đã làm cho tài liệu có nội dung phù hợp hơn, tránh được nhiều sai sót. Xin chân thành cám ơn về sự đóng góp quí giá trên. Người viết căn cứ sự giúp đỡ, chỉ dẫn trên đã cố gắng hoàn chỉnh tài liệu kịp ra mắt bạn đọc, tuy nhiên vì khả năng có hạn chắc rằng trong tài liệu vẫn còn những sai sót khó tránh. Rất mong bạn đọc gần, xa góp thêm ý kiến nhằm làm cho tài liệu ngày càng hoàn thiện. Thư, bài góp ý, xây dưng xin gửi về phân hiệu Đại học Hàng Hải, thành phố Hồ Chí Minh. Thành phố Hồ Chí Minh tháng 12 năm 2000. 5
  7. Chương 1 PHƯƠNG PHÁP TÍNH VÀ TỰ ĐỘNG HÓA TÍNH TOÁN, THIẾT KẾ TÀU Trong chương này giới thiệu những phương pháp tính thông dụng dùng xử lý những vấn đề thường gặp trong tính tính nổi tàu thủy, tính ổn định tàu, thiết kế máy đẩy tàu, thiết kế tối ưu tàu thủy. 1.1 NỘI SUY LAGRANGE Đa thức nội suy Lagrange được viết dưới dạng1: f(x) = pn(x) + Rn(x), (1.1) hoặc dạng đầy đủ: n ⎡ n ⎤ f (n+ 1) ()ξ f()()()() x= Li x f xi + Π x − xi ,a<ξ < b (1.2) ∑ ⎢i=0 ⎥ i=0 ⎣ ⎦ (n + 1)! n ⎛ x− x ⎞ trong đó L() x = ⎜ j ⎟ (1.3) i ∏⎜ ⎟ j=0 ⎝ xi− x j ⎠ j≠ i n Đa thức pn ()()() x= ∑ Li x f xi mang tên gọi đa thức Lagrange, còn vế sau của phía phải công thức i=0 gọi hàm sai số. Đa thức pn(x) mặt khác được hiểu là đa thức bậc n, có dạng: pn(x) = a0(x – x1) (x – x2) (x – xn) + + a1(x – x0) (x – x2) (x – xn) + + a2(x – x0) (x – x1) (x – xn) + + ai(x – x0) (x – x1) (x – xi-1) (x – xI+1) (x – xn) an(x – x0) (x – x1) (x – xn-2)(x – xn-1 ) (1.4) Các hệ số a0, a1, a2, tính từ quan hệ: pn(xi) = f(xi) = yi ; i = 0, 1, 2, (1.5) Lần lượt thay x = x0, x = x1, vào công thức cuối có thể xác định công thức tính các hệ số. Ví dụ, từ pn(x0) = y0 = a0(x0 – x1) (x0 – x2) (x0 – xn) sẽ nhận được: f() x0 a0 = (x0− x 1 )( x 0 − x 2 ) ( x 0 − xn ) tương tự vậy có thể viết: 1 R.W. Hamming, “Numerical Methods for Scientists and Engineers”, McGraw-Hill, N.Y, 1962, F.B. Hildebrand, “Introduction to Numerical Analysis”, McGraw-Hill, N.Y., 1956. 6
  8. f() x1 a1 = (x1− x 0 )( x 1 − x 2 ) ( x 1 − xn ) f() xn an = (xn − x0 )( xn − x1 ) ( xn − x n−1 ) Hệ số thứ i mang dạng chung: f() xi ai = (1.6) (xi − x0 )( xi − x1 ) ( xi − x i−1 )( x i − x i+1 ) ( xi− x n ) Thay các biểu thức vừa xác định vào vị trí a0, a1, , an sẽ nhận được công thức nội suy hay còn gọi đa thức Lagrange: (x− x )( x − x ) ( x− x ) (x− x )( x− x ) ( x− x ) p() x = 1 2 n y + 0 2 n y + n (x− x )( x − x ) ( x − x ) 0 (x− x )( x − x ) ( x − x ) 1 0 1 0 2 0 n 1 0 1 2 1 n (1.7) (x− x0 )( x − x1 ) ( x − xn−1 ) + + yn (xn − x0 )( xn − x1 ) ( xn − x n−1 ) n n ⎛ x− x ⎞ hoặc dưới dạng gọn hơn như đã trình bày p()()() x= L x f x , với L() x = ⎜ j ⎟ . n ∑ i i i ∏⎜ ⎟ i=0 j=0 ⎝ xi− x j ⎠ j≠ i Những trường hợp riêng lẻ của hàm nội suy Lagrange như sau. với n =1: ()x− x1 ()x− x0 p1 () x = y0 + y1 (1.8) ()x0− x 1 ()x1− x 0 với n = 2: (x− x1 )( x − x2 ) (x− x0 )( x − x2 ) (x− x0 )( x − x1 ) p2 () x = y0 + y1 + y2 (1.9) (x0− x 1 )( x 0 − x 2 ) (x1− x 0 )( x 1 − x 2 ) (x2− x 0 )( x 2 − x 1 ) Hàm p1(x) là đoạn thẳng qua hai điểm (x0, y0) (x1, y1 ), có tên gọi công thức nội suy tuyến tính. Hàm thứ hai là đường parabol bậc hai qua ba điểm cho trước, gọi là nội suy bậc hai. Chương trình hóa phương pháp nội suy Lagrange được thể hiện bằng thuật toán Neville và minh họa tại hàm bằng ngôn ngữ C sau. #include (math.h> void Lagrange(xa, ya, n, x, y, dy) float xn[], ya[], x,*y, *dy; int n; { int i, m, ns=1; float den, dif, dift, h0, hp, w; float *c, *d, *vector(); dif = fabs( x-xa[q]); c = vector(1,n); d = vector(1,n); for (i=1; i<=n; i++) { if ( ( dift = fabs(x - xa[1] )) < dif) { ns =i; 7
  9. dif = dift; } c[i] = ya[i]; d[i] = ya[i]; } *y = ya[ns ]; for ( m=1; m<n; m++) { for (i=1; i <= n-m; i++) { h0= xa[i] -x; hp = xa[i+m] - x; w = c[i+1] - d[i]; if ((den=h0-hp) == 0.0) nerror("Error here !"); den = w/den; d[i] = hp*den; c[i] = h0*den; } *y += (*dy=(2*ns < (n-m) ? c[ns+1]: d[ns ] )); } free_vector(d,1,n); free_vector(c,1,n); } 1.2 TÍCH PHÂN MỘT LỚP Giả sử cần thực hiện tích phân hàm f(x) từ a đến b, có thể tiến hành tích từng phân đoạn và sau đó tổng hợp kết quả. Đoạn [a, b] được chia làm n phân đoạn bằng nhau, giới hạn bằng các nút, trong đó nút đầu tiên x0 = a, tiếp đó x1, x2, , xn-1, xn= b. Giá trị hàm f(x) được xác định cho tất cả các nút. b −1 Nếu ký hiệu: h = =x − x , mang tên gọi bước, có thể thay biến x bằng biến mới như sau: n k+1 k x− x p = k . h Hàm f(x) giờ có thể viết: p( p − 1) 2 f(x) = f(x0 + ph) = f(x0) + p.Δf(x0) + Δ f() x + + R(x0 + p.h) (1.10) !2 0 trong đó hàm sai số được tính theo công thức: n+1 (n+1) R(x0 + ph) = h p(p-1) (p-n)/ (n+1)! * [f (ξ) / 2! ], ξ ∈(x0, xn). (1.11) Công thức tính tích phân nêu trên có thể viết thành: (1) Trong đa thức cuối chỉ giữ lại một thành phần đầu tiên, kết quả nhận được công thức tính tích phân theo nguyên tắc hình chữ nhật. b n−1 xk +1 k= n −1 f() x dx≅ f() x dx= hf( x )= h f ( x ) + + f ( x ) ∫ ∑ ∫ k ∑ k [0 n−1 ] k=0 k =0 a xk (2) Nếu giữ lại hai thành phần đầu của biểu thức, kết quả sẽ nhận được công thức tính tích phân theo luật hình thang: 8
  10. b n−1 xk +1 k= n −1 1 n−1 f()() x+ f x f() x dx≅ f() x dx= h[ f ( x )+ p Δ f ( x )] du = h k k +1 = ∫ ∑ ∫ k ∑ ∫ k ∑ k =0 k =0 k =0 2 a xk 0 = h[] f( x0 ) / 2+ f ( x1 ) + + f ( xn−1 ) + f ( x2 ) / 2 (3) Công thức Simpson được tính theo luật trên đây khi giữ lại ba thành phần đầu tiên của chuỗi. b 2n− 2 xk +2 2n− 2 2 p( p − 1) f() x dx ≅ f() x dx= h[()() f x+ p Δ f x + Δ2 f( x )] du = ∫ ∑ ∫ k ∑ ∫ k k k k=0,2, k=0,2, 2! a xk 0 h 2n− 2 h = ∑[f ( xk )+ 4 f ( xk+1 ) + f ( xk+2 )] =[]f( x0 ) + 4 f ( x1 ) + 2 f ( x2 ) + + 4 f ( x2n− 1 ) + f ( x2n ) 3 k=0,2, 3 Hàm bằng ngôn ngữ C, thực hiện tích phân theo phương pháp hình thang được trình bày tiếp dưới đây. #define FUNC(x) ((*func) (x) ) float trapezd( func, a, b, n) float a, b; float (*func) (); int n; { float x, tnm, sum, del; static float s; static int it; int j; if (n == 1 ) { it=1; return (s=0.5* (b-a) *(FUNC(a) +FUNC(b) )); } else { tnm = it; del = (b-a)/tnm; x = a + 0.5*del; for (sum=0.0, j=1; j<it; j++, x+=del) sum += FUNC(x); it *= 2; s=0.5* (s + (b-a)*sum / tnm); return s; } } Trường hợp khoảng chia theo trục Ox không bằng nhau, công thức tính theo đề nghị của Milne 2 được sử dụng. Tích phân hàm f(x) trong phân đoạn xác định bằng ba nút được tính theo đề nghị của Milne: x3 f()()()() x dx= a f x + a f x + a f x (1.12) ∫ i 1i 1 2i 2 3i 3 x1 2 Nếu gán fi(x) các giá trị 1, x và x , các ẩn số ai được xác định theo công thức: 2 Milne, W.E. “Numerical Calculus”, Princeton, 1949 9
  11. ⎡1 1 x3− x 1 ⎤ a3=() x 3 − x 1 ⎢ − ⎥; ⎣2 6 x3− x 2 ⎦ 3 1 ()x3− x 1 a2 = ; (1.13) 6 (x2− x 1 )( x 3 − x 2 ) a1= x 3 − x 1 − a 2 − a 3 Hàm viết bằng ngôn ngữ C xử lý tích phân giới hạn biến thiên được thể hiện như sau. x integrationv = ∫ f() x dx 0 void integrationv(int N, float *X, float *f, float *Ans) { register i,ii,n2 ; float d1,d2,d3,a1,a2,a3 ; n2 = N / 2; Ans[0]=0.0; for (i = 1; i <= n2; i++) { ii = 2*i-1; d1=X[ii]-X[ii-1]; Ans[ii]=Ans[ii-1]+d1/2.0*(f[ii]+f[ii-1]); if (ii != N-1) { d2=X[ii+1]-X[ii-1]; d3= d2/d1; a2= d3/6.0*sqr(d2)/(X[ii+1]-X[ii]); a3= d2/2.0-a2/d3; Ans[ii+1]=Ans[ii-1]+(d2-a2-a3)* f[ii-1]+a2*f[ii]+a3*f[ii+1]; } } } Cách làm trên đây được áp dụng tính gần đúng tích phân, trong đó các bước chia trên trục không bắt buộc đều nhau, số phân đoạn lớn hơn 1, được thể hiện tại hàm sau đây: float integ(int N, float *X, float *f ) { int M, Lb, k ; float Sum ; M=N-1; Sum=(X[1]-X[0])/6.0*(f[0]*((X[1]-X[2])/(X[0]-X[2])+2.0)+ f[1]*((X[0]-X[2])/(X[1]-X[2])+2.0)-f[2]*sqr(X[1]-X[0])/ ( (X[0]-X[2])*(X[1]-X[2]) ) ); Lb=1; if ( N != 3 ) { Sum=Sum+(X[2]-X[1])/6.0*(f[1]*((X[2]-X[3])/(X[1]-X[3])+2.0)+ f[2]*((X[1]-X[3])/(X[2]-X[3])+2.0)-f[3]*sqr(X[2]-X[1])/ ((X[1]-X[3])*(X[2]-X[3]))); Lb=2; } for (k = Lb ; k < M; k++) Sum=Sum+(X[k+1]-X[k])/6.0*(f[k]*((X[k+1]-X[k-1])/(X[k]-X[k-1])+2.0) 10
  12. + f[k+1]*((X[k]-X[k-1])/(X[k+1]-X[k-1])+2.0)-f[k-1]* sqr(X[k+1]-X[k])/((X[k]-X[k-1])*(X[k+1]-X[k-1]))); return Sum; } 1.3 ĐA THỨC LEGENDRE Đa thức Legendre trực giao trong đoạn [-1, 1], được hiểu như sau: +1 P( x ) P ( x ) dx= 0, n ≠ m (1.14) ∫ n m −1 +1 P( x ) P ( x ) dx= c ( n ) ≠ 0 ∫ n n −1 Một số ít đa thức Pn (x) được viết dưới đây: P0(x) =1, P1(x) =x, 2 P2(x) = ½ (3x - 1), 3 P3(x) = ½ (5x – 3x), P(x) = 1/8 (35x4 – 30x2 + 3) Quan hệ hồi qui của đa thức Legendre có dạng: 2n − 1 n −1 P() x = xP() x − P() x . (1.15) n n n−1 n n−2 Áp dụng tính trực giao của đa thức Pn(x) để xác định các tọa độ xj nhằm giảm thiểu sai số khi tính. Công thức tính tích phân Gauss-Legendre có dạng: b b b f()()() x dx= p x dx+ R x dx (1.16) ∫∫n ∫ n a a a Trong đó hàm f(x) có thể viết: f(x) = pn(x) + Rn(x). n ⎡ n ⎤ f (n+ 1) ()ξ f()()()() x= Li x f xi + Π x − xi ,a<ξ < b (1.17) ∑ ⎢i=0 ⎥ i=0 ⎣ ⎦ (n + 1)! n ⎛ x− x ⎞ trong đó L() x = ⎜ j ⎟ i ∏⎜ ⎟ j=0 ⎝ xi− x j ⎠ j≠ i Nếu sử dụng biến z theo cách sau đây: a+ b b− a z = + x (1.18) 2 2 n ⎛ z− z ⎞ hàm f(x) có thể thay bằng F(z), hàm Lagrangre L (x) có thể thay bằng L() z = ⎜ j ⎟ còn i i ∏⎜ ⎟ j=0 ⎝ zi− z j ⎠ j≠ i giới hạn tích phân trở thành –1 và +1, biến ξ nằm trong giới hạn -1 < ξ < +1. Tích phân (1.19 ) trở thành: 11
  13. +1 +1 n +1⎡ n ⎤ F() z dz= L()()()() z F z dz +z − z q z dz (1.19) ∫ ∫∑ i i ∫ ⎢∏ i ⎥ n −1 −1 i=0 −1⎣ i=0 ⎦ Nếu bỏ vế sau bên phía phải của công thức có thể viết: +1 n +1 n F()()() z dz≅ F z L z dz= W F() z (1.20) ∫ ∑ i∫ i ∑ i i −1 i=1 −1 i=1 Từ công thức cuối có thể xác định: +1 +1 (z− z0 ) ( z − zi−1 )( z − zi+1 ) ( z − zn ) dx w= L() z dz = ∫−1 (1.21) i∫ i −1 (z− z0 ) ( z − zi−1 )( z − zi+1 ) ( z − zn ) +1⎡ n ⎤ Vế thứ hai của phía phải ()()z− z q z dz chính là sai số của phép tích phân số đang được xét. ∫ ⎢∏ i ⎥ n −1⎣ i=0 ⎦ Trong giai đoạn này cần thiết xác định vị trí của zi trong đa thức Legendre nhằm làm cho vế này trượt tiêu. Tại đây sử dụng tính trực giao của hàm Legendre để đạt điều momg muốn. Khai triển n hai đa thức qn(z) và ∏()z− zi dưới dạng sau đây: i=0 n n+1 ∏(z− zi ) = b0 P 0 ( z ) + b1 P 1 ( z ) + = ∑bi P i ( z ) (1.22) i=0 i=0 và n q(z) = c0P0(z) +c1P1(z) + = ∑ci P i () z (1.23) i=0 Thay hai đa thức cuối vào biểu thức thuộc vế hai phía phải công thức (1.19) có thể viết: +1⎡ n n n ⎤ b c P()() z P z+ b c P()() z P z dz (1.24) ∫ ⎢∑∑ i j i j n+1 ∑ i i n+1 ⎥ −1⎣ i==00j i=0 ⎦ Từ tính trực giao của đa thức Legendre có thể viết: +1 b c P( z ) P ( z ) dz= 0; i ≠ j (1.25) i j∫ i j −1 Từ đó, sai số phép tích phân được viết là: +1⎡ n ⎤ +1 n n +1 ()()z− z q z dz= b c P() z2 dz= b c P() z2 dz (1.26) ∫ ⎢∏ i ⎥ n ∫∑ i i[] i ∑ i i ∫ []i −1⎣ i=0 ⎦ −1 i=0 i=0 −1 Từ phương trình cuối có thể xác định các nút tính toán. Các nút tìm bằng cách này gọi là các điểm zero của đa thức Legendre. b Giới hạn a, b trong công thức chung f() z dz , có quan hệ với giới hạn chuẩn của hàm f(x) như đã ∫a a+ b b− a trình bày z = + x . 2 2 Công thức tính tích phân theo Gauss-Legendre có dạng: 12
  14. b n f()()() x dx≅ b − a W f z (1.27) ∫ ∑ i i a i=1 Một số giá trị của nút và hệ số tính cho công thức trên như sau: Bảng 1.1 n=1 x1= 0,0 w1 = 2,0 n=2 x1,2 = ±0,577350 a1,2 = 1,0 n=3 x3,1 = ±0,774597 a1,3 = 0,555556 x2 = 0,0 a2 = 0,888889 n=4 x3,1 = ±0,861136 a1,4 = 0,347855 x4,2 = 0,339981 a2,3 = 0,652145 Ví dụ tính theo công thức Gauss-Legendre. 2 dx Tính giá trị tích phân sau ∫1 x 2 dx Thay thế biến z = 2x – 3 và hàm f(x) thành F(z) = 2/ (z+3). Tích phân có giá trị bằng tích ∫1 x +1 2 phân dz . Kết quả tính, với n = 4 sẽ có dạng: ∫−1 z + 3 Bảng 1.2 i zi wi F(zi) wF(zi) 0 0,0 0,5688889 0,333333 0,189629 1 0,5384693 0,4786287 0,282608 0,135264 2 -0,5384693 0,4786287 0,406251 0,1944435 3 0,9061798 0,2369269 0,256004 0,0606544 4 -0,9061798 0,2369269 0,4775959 0,1131553 Cộng: 0,6931471 Hàm viết bằng ngôn ngữ C cho trường hợp n =10 được ghi lại dưới đây. float qgaus (func, a, b) /* Gauss-Legendre formula */ float a,b; float (*func) (); { int j; float xr, xm, dx, s; static float x[] = {0.0, 0.1488743389, 0.4333953941, 0.679409568, 0.8650633666, 0.9739065285}; static float w[] = { 0.0, 0.2955242247, 0.269266719, 0.2190863625, 0.149451349, 0.0666713443}; xm=0.5* (b+a); xr=0.5*(b-a); s=0.0; for (j=1; j<=5; j++) { dx=xr*x[j]; s += w[j]* (( *func)(xm+dx) + (*func) (xm-dx) ); 13
  15. } return s *= xr; } 1.4 ĐA THỨC TCHEBYSHEV Đa thức Tchebyshev, ký hiệu từ ký tự đầu tên nhà bác học Tn(x), trực giao trong [-1, 1] cùng hàm trọng lượng w(x) = 1/ √ (1 – x2). +1 1 Tn( x ) T m ( x ) dx= 0, n ≠ m (1.28) ∫ 2 −1 1− x +1 1 Tn( x ) T n ( x ) dx= c ( n ) ≠ 0 (1.29) ∫ 2 −1 1− x Những đa thức đầu có dạng sau. T0(x) =1, T1(x) =x, 2 T2(x) = 2x - 1, 3 T3(x) = 4x – 3x. Quan hệ hồi qui của đa thức Tchebyshev có dạng: T(x) = 2xTn-1(x) – Tn-2(x) (1.30) Công thức Gauss – Tchebyshev. Tích phân hàm f(x), giới hạn a, b theo công thức Gauss – Tchebyshev sẽ là: b ()b− a + ⎛ a + b b− a ⎞ ∫ f() x dx = ∫ f ⎜ + z⎟ dz (1.31) a 2−1 ⎝ 2 2 ⎠ Hàm trọng lượng được trình bày trên, cho phép viết tiếp: b ()b− a + 1 f() x dx = F() z dz , trong đó ∫ ∫ 2 a 2 −1 1− z ⎛ a+ b b− a ⎞ F( z )= 1 − z2 f ⎜ + z ⎟ (1.32) ⎝ 2 2 ⎠ và b π n f() x dx ≅(b − a ) 1 − z2 f ( x ) (1.33) ∫ ∑ i a 2n i=1 a+ b b− a x = + z i 2 2 i ⎛ []2(i − 1) + 1π ⎞ zi = cos⎜ ⎟ ⎝ 2n ⎠ Các nút trong công thức Tchebyshev được trình bày tại bảng 1.3. 14
  16. Bảng 1.3 n Wi X 3 2/3 X1 = -X3 = 0,707107 X2 = 0,0 4 ½ X1 = -X4 = 0,794654 X2 = -X3 = 0,187592 5 2/5 X1 = -X5 = 0,832498 X2 = -X4 = 0,374541 X3 = 0,0 1.5 TÌM NGHIỆM PHƯƠNG TRÌNH BẰNG PHƯƠNG PHÁP CHIA ĐÔI ĐOẠN CÓ NGHIỆM Phương pháp tìm nghiệm phương trình y = f(x) giản đơn mà nhanh, được dùng tại đây mang tên gọi bisection ( lưỡng đoạn). Giả sử phương trình f(x) = 0, trong đó f(x) liên tục trong đoạn [a, b] và f(a).f(b) #define Max 40 15
  17. float rtbis( func, x1, x2, xacc) float x1, x2, xacc; float (*func) (); { int j; float dx, f, fmid, xmid, rtb; void nerror(); f=(*func) (x2); if (f*fmid >= 0.0) nerror("Bisection Method); rtb = f < 0.0 ? (dx=x2-x1, x2); for (j=1; j <= Max; j++) { fmid=(*func) (xmid=rts+(dx *= 0.5) ); if (fmid <= 0.0) rtb=xmid; if (fabs(dx) <xacc || fmid == 0.0) return rtb; } nerror("Too many bisections "); } 1.6 PHƯƠNG PHÁP TỔNG NHỎ NHẤT CÁC BÌNH PHƯƠNG Đây là phương pháp kinh điển, ra đời từ thế kỷ thứ XIX áp dụng cho trường hợp hàm hóa, biểu diễn kết quả thực nghiệm, đo đạc thực tế. Giả sử rằng từ dữ liệu thực tế có thể xác lập quan hệ giữa hai biến x và y như sau. Ứng với mỗi giá trị cửa xj, j =0,1, 2, , m có thể nhận biết giá trị thực của y là yj = f(xj). Nếu giờ đây cần phải xác lập một hàm xấp xỉ biểu diễn quan hệ đó dưới dạng : p(x) = c0f 0(x) + c 1f 1(x) + c 2f 2(x) + + cnf n(x) (a) trong đó n < m. Hàm fk(x), k =0, 1, 2, , n liên tục chứa biến x, còn ck là ẩn số phải xác định. Trường hợp k riêng, song thường gặp, fk(x) = x , hàm p(x) cần xác định chính là đa thức bậc n. Theo cách tổ chức này, sai số xác định tại mỗi vị trí xj, j =0, 1,, , m giữa hàm p(x) và f(xj) được tính như sau: p(x) – f(xj) = δ j j = 0, 1, 2, , m (b) Phương pháp tính đòi hỏi tổng bình phương các sai số phải nh ỏ nhất : m m 2 2 ∑δ j = ∑[]p( x ) − f ( x j )→ min (c) j=0 j=0 Điều kiện này cho phép viết : m m ∂ 2 ∂δ j ∑δj = ∑ δ j =0;k = 0,1,2, ,m (d) ∂ck j=0 j=0 ∂ck ∂δ j với = f() x j ∂ck Dưới dạng đầy đủ công thức cuối trở thành: m ⎡ n ⎤ ∑∑⎢ []ci f i( x j )− y j ⎥ fk ( x j )= 0; k = 0,1,2, ,m (e) j==00⎣ i ⎦ hoặc có thể viết dưới dạng, thuận lợi hơn khi lập trình: 16
  18. n ∑ aki ci= b k ; k = 0,1,2, , m ( ) f i=0 m aki = ∑ fi( x j ) f k ( x j ); j=0 với m bk = ∑ yj f k() x j j=0 x 2x Ví dụ: Xác lập hàm y = c1 e + c2e qua ba điểm (0,1) (1, -2) (2, -40) trên cơ sở phương pháp tổng nhỏ nhất các bình phương sai số. 0 0 c1 e + c2 e - 1 = δ1 1 2 c1 e + c2 e - 1 = δ2 2 4 c1 e + c2 e - 1 = δ3 Theo cách làm nêu trên, có thể viết: ∂δ ∂δ ∂δ δ 1 + δ 2 + δ 3 = 0 1 ∂c2 ∂ c3 ∂ c 1 1 1 ∂δ 1 ∂δ 2 ∂δ 3 δ 1 + δ 2 + δ 3 = 0 ∂c2 ∂ c2 ∂ c2 ∂δ∂ δ ∂δ ∂δ∂ δ ∂δ Thay các giá trị 1 = e0 ;;2 = e1 3 = e 2 ; 1 = e0 ;;2 = e2 3 = e4 vào hệ ∂c1 ∂c1 ∂c1 ∂c2 ∂c2 ∂c2 phương trình cuối có thể xác định : c1 ≈ 2; c2 ≈ -1. Nếu sử dụng các ký hiệu vector và ma trận vào bài toán này, hệ phương trình nêu trên được viết lại gọn hơn. Các hệ số c sắp xếp tại vector {c} gồm n+1 thành phần, giá trị yj, j = 0, 1, 2, , m xếp vào {y}, ma trận [F] tập họp m+1 dòng, ứng với các gía trị của fk(xj) . ⎧c0 ⎫ ⎧ y0 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪c1 ⎪ ⎪ y1 ⎪ {}c = ⎨ ⎬ ;{y } = ⎨ ⎬ ; ⎪ • ⎪ ⎪ • ⎪ ⎪c ⎪ ⎪y ⎪ ⎩ n ⎭(N+ 1) x 1 ⎩ m ⎭(M+ 1) x 1 ⎡ f0()()() x 0 f 1 x 0 • fn x m ⎤ ⎢ f()()() x f x• f x ⎥ ⎢ 0 1 1 1 n m ⎥ ⎢ •••• ⎥ ⎢ ⎥ ⎢ •••• ⎥ ⎢ f()()() x f x• f x ⎥ ⎣ 0 m 1 m n m ⎦ (M+ 1) x ( N + 1) Vector sai số xác định bằng biểu thức: [F]{c} – {y} = {δ}, Và công thức (d) trở thành: 17
  19. ∂ 2[] [F ]{ c }− { y } [ F ]{ c } = 0 ∂{}c hoặc : [][F ]{ c }− { y }[ F ] = 0 Vì rằng [F] ≠ {0}, biểu thức trong ngoặc vuông phải bằng 0. Bài toán theo phương pháp tổng nhỏ nhất các bình phương đưa về dạng: [F](m1xn1).{c}(n1x1) = {y}(m1x1), trong đó m1 = m + 1; n1 = n + 1. (g) Ma trận [F] gồm m+1 dòng, n+1 cột, với n < m. Vector {c} chỉ có n+1 d ò ng, {y } gồm m+1 dòng. Xử lý phương trình (g) theo nhiều cách khác nhau. Dưới đây trình bày hai cách thông dụng, dễ dùng. Cách thứ nhất đưa bài toán (g) về dạng bài toán kinh điển như trình bày tại (f). Nhân hai vế của phương trình (g) với ma trận chuyển vị của [F] là [F]T : [F]T[F]{c} =[F]T{y} Từ đó có thể viết : [A].{c} = {b}, (h) với [A] = [F]T[F] và {b} = [F]T{y}. (i) Các hệ số {c} được xác định sau khi giải hệ phương trình đại số tuyến tính. Ví dụ : áp dụng cách làm này xử lý bài toán vừa nêu. ⎡e0 e 0 ⎤ ⎧ 1 ⎫ ⎢ 1 2 ⎥⎧c1 ⎫ ⎪ ⎪ e e ⎨ ⎬ =⎨ − 2 ⎬ ⎢ ⎥ c ⎢ 2 4 ⎥⎩ 2 ⎭ ⎪ ⎪ ⎣e e ⎦ ⎩− 40⎭ ⎡e0 e 0 ⎤ 1 0 1 2 0 1 2 ⎧ ⎫ ⎡e e e ⎤⎢ 1 2 ⎥⎧c1 ⎫ ⎡e e e ⎤⎪ ⎪ ⎢ ⎥ e e ⎨ ⎬ = ⎢ ⎥⎨ − 2 ⎬ e0 e 2 e 4 ⎢ ⎥ c e0 e 2 e 4 ⎣ ⎦⎢ 2 4 ⎥⎩ 2 ⎭ ⎣ ⎦⎪ ⎪ ⎣e e ⎦ ⎩− 40⎭ ⎡ 63 424 ⎤⎧c1 ⎫ ⎧ − 296 ⎫ ⎢ ⎥⎨ ⎬ = ⎨ ⎬ ⎣424 3036⎦⎩c2 ⎭ ⎩− 2198⎭ Từ hệ phương trình cuối có thể xác định : c1 ≈ 2; c2 ≈ -1. Cách làm thứ hai là tiến hành giải hệ phương trình đại số gồm m+1 phương trình, chứa chỉ n +1 ẩn, trong đó n < m, bằng phương pháp xử lý ma trận suy biến. Phương pháp mang tên gọi bằng tiếng Anh: Singular Value Decomposition – SVD. Một trong những thuật toán hay nhất xử lý hệ phương trình đại số tuyến tính [A] (MxN){b{(Nx1)={y}(Mx1) được trình bày trong sổ tay toán tính của Wilkinson. Thủ tục tính được nhóm nhà toán học dưới sự chỉ dẫn của Wilkinson viết ra trong những năm sáu mươi bằng ngôn ngữ Algol đã dùng có hiệu quả trên những dàn máy tính. Thủ tục trên người viết tài liệu này đã “dịch” sang ngôn ngữ FORTRAN và sau đó “chuyển “ sang ngôn ngữ C, được giới thiệu tiếp theo, giúp bạn đọc xử lý những bài toán thực tế. Trong chương trình, người viết đang hạn chế n <10. Trường hợp 18
  20. người dùng cần tăng số hệ số ck đến số lớn hơn 9 ( tính từ 0), cần thay thế JMAX bằng số hệ số cao nhất. Hàm LeastSquares được dùng xác định các hệ số thủy động lực chân vịt tàu, xác định các hệ số đường cong sức cản vỏ tàu trong phần tự động hoá thiết kế chân vịt. Hàm này là phương tiện chính cho các phép hồi qui dùng trong phần xác định sức cản tàu. #include #include #include #include #define JMAX 10 void LeastSquares( int N, int M, double A[][JMAX], double *b, double *y ) { register int N1, M1, MM1, L, i, LL, L1, j; double s, ss, s2, d, pp; MM1 = M -1; M1 = M+1; for ( L =0; L < M; L++) { ss =0.0; for ( i = L; i < N; i++) ss += A[i][L] * A[i][L]; s2 = ss; s = sqrt(s2); if ( A[L][L] < 0.0 ) s = -s; d = s2 + s*A[L][L]; A[L][L] += s; if ( L != MM1) { L1 = L + 1; for ( j =L1; j < M; j++) { pp = 0.0; for ( i =L; i < N; i++) pp += A[i][L] * A[i][j]; A[N][j] = pp/ d; } for ( j = L1; j < M; j++) for ( i =L; i < N; i++) A[i][j] -= A[i][L] * A[N][j]; } A[N][L] = -s; } for ( i=0; i <N; i++) A[i][M] = y[i]; for ( L =0; L < M; L++) { pp=0.0; for ( i=L; i < N; i++) pp += A[i][L] * A[i][M]; d = pp / ( -A[L][L] * A[N][L] ); for ( i= L; i < N; i++) A[i][M] -= d * A[i][L]; } b[MM1] = A[MM1][M] / A[N][MM1]; if ( N != 1) { for ( LL =0; LL < MM1; LL++) { L = M - LL -2; L1 = L+1; pp = A[L][M]; for ( i=L1; i < M; i++) pp -= A[L][i] * b[i]; b[L] = pp / A[N][L]; } } 19
  21. s =0.0; for ( i=M; i < N; i++) { ss += A[i][M] * A[i][M]; A[i][M] = 0.0; } for ( LL =0; LL < M; LL++) { L = M-LL-1; pp =0.0; for ( i = L; i < N; i++) pp += A[i][L] * A[i][M]; d = pp / ( -A[L][L] * A[N][L] ); for ( i=0; i <N; i++) A[i][M] -= d * A[i][L]; } return; } main() { int m, n, Npoints,i, j, jp1,Min, Max; double a[200][JMAX], b[JMAX]; double xi,p, res, ss; int order, coeff[10] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.0 }; double x[10] ={ 0.23, 1.01, 2.29, 2.87, 4.15, 5.36, 5.51, 6.36, 6.84, 7.0}; double y[10]= {5.64, 7.83, 17.04, 21.38, 24.56, 16.21, 14.57, 0.78, -7.64, -12.52 }; printf("Enter n, npoints\n"); sca nf ("%d %d", &n, &Npoints); for ( m =0; m < Npoints; m++) { a[m][0] = 1.0; xi = x[m]; p =1.0; for ( j = 0; j < n; j++) { jp1 =j + 1; p *= xi; a[m][jp1] = p; } } Least Squares( Npoints, n+1, a, b, y ); res =0.0; printf("\n i x[i] y[i] observation\n"); for ( m=0; m< Npoints; m++) { p =1.0; xi = x[m]; ss =0.0; for ( j =0; j < n + 1; j++) { ss += b[j] * p; p *= xi; } printf(" %d %12.6lf %15.8lf %15.8lf\n", m, x[m], y[m], ss); ss -= y[m]; ss = ss * ss; res += ss; } printf("\n\nThe order is %d with res = %10.6lf\n\n", n, res); for ( i =0; i <= n; i++) printf("\tC[%d] = %15.8lf\n", i, b[i]); } 20
  22. 1.7 QUI HOẠCH TUYẾN TÍNH Trong quản lý sản xuất thường gặp những vấn đề không đơn giản cần được cân nhắc, tính toán thiệt hơn trên cơ sở nguồn nhân lực và tài nguyên đang sở hữu, trước khi đưa ra những quyết định. Tài nguyên trong mỗi xí nghiệp có khác nhau, tuy nhiên những tài nguyên cấp thiết nhất mà các xí nghiệp sản xuất phải có gồm máy móc, thiết bị, lao động, tiền bạc, nguyên vật liệu và nhà xưởng, văn phòng. Qui hoạch tuyến tính trong tài liệu này là phương pháp tính, tận dụng nguồn nhân lực và tài nguyên, chọn cách làm hữu hiệu nhất cho công việc. Qui hoạch tuyến tính đóng vai trò phương pháp tính toán để giúp cho người điều hành sản xuất tận dụng đến mức hợp lý nhất tài nguyên hiện hữu, đưa ra những quyết định hợp thời cho sản xuất nhằm đưa lại hiệu quả lớn nhất. Theo nghĩa này qui hoạch tuyến tính được hiểu là cách tìm maximum hoặc minimum của bài toán đang đặt ra cho sản xuất, trong không gian hạn chế và thời gian cũng hạn định. Dạng chung của bài toán được phát biểu theo cách sau. Tìm maximum (hoặc minimum) Z = c1X1 + c2X2 + + cnXn (1.34) Với điều kiện: A11X1 + a12X2 + + a1nXn ≤ b1 A21X1 + a22X2 + + a2nXn ≤ b1 am1X1 + am2X2 + + amnXn ≤ b1 (1.35) và hạn chế: X1 ≥ 0; X2 ≥ 0 ; Xn ≥ 0. (1.36) Dưới dạng vector bài toán trên được thể hiện: Hàm mục tiêu: cTX → maximum (1.37) Hạn chế: AX = b; (1.38) X ≥ 0 (1.39) Ví dụ đơn giản sau đây minh hoạ cách phát biểu trên lấy từ xí nghiệp sản xuất đồ mộc. Xí nghiệp cần sản xuất lượng bàn và ghế nhằm thu lợi nhuận cao nhất. Thực tế sản xuất cho biết rằng mỗi bàn đòi 4 giờ công mộc, 2 giờ công đánh véc ni, sơn. Công sản xuất ghế là 3 giờ mộc, 1 giờ đánh vecni và sơn phết. Xí nghiệp giành 240 công làm mộc và 100 công sơn phết cho công việc này. Tiền lời do bán bàn 7.000 Đ, lời từ bán ghế 5.000 Đ. Theo cách viết trên đây, X1 dùng chỉ số bàn sẽ sản xuất, X2 – số ghế. Giá bán cả bàn và ghế 3 3 tính bằng VNĐ, thời điểm 1998 là 7.10 .X1 + 5.10 X2. Bài toán qui hoạch theo mô hình vừa xây, tính bằng đơn vị 103 VNĐ: 7. X1 + 5. X2 → maximum. (a) 21
  23. 4XX1+ 3 2 ≤ 240⎫ ⎬ ( ) b 2XX1+ 1 2 ≤ 100⎭ và: X1 ≥ 0; X2 ≥ . 0 ( ) c Cách giải kinh điển cho bài toán này có thể kể: • phương pháp đồ thị, • phương pháp iso-Profit (đẳng lợi nhuận) • phương pháp chọn điểm tối ưu. Trong tài liệu này chỉ trình bày phương pháp thích hợp cho việc tự động hóa, gọi là phương pháp “giản đơn “ (Simplex Method). Thủ tục thực hiện trải qua các bước. Bước đầu tiên xử lý các điều kiện hạn chế đã đặt ra trong bài toán, chuyển hoá các bất đẳng thức thành các đẳng thức, bằng cách gán thêm các biến lỏng lẻo (mượn cách gọi của người Mỹ: slack variables). 4XXS1+ 3 2 + 2 = 240⎫ ⎬ ( ) ’ b 2XXS1+ 1 2 + 1 = 100 ⎭ Viết đầy đủ hơn, hệ phương trình trên sẽ là: 4XXSS1 + 32 + 0 1 + 2 = 240⎫ ⎬ ) ( ’ ’ b 2XXSS1+ 1 2 + 1 1 + 0 2 = 100⎭ Hàm mục tiêu trở thành: . 7 X1 + 5 m X2 +0S( 1 +0S2 → u ) m ’ i a x a m Tiếp đó áp dụng phép thử và giải (b’’), bắt đầu từ X1 = X2 = 0 và lợi nhuận bằng 0. Trong trường hợp này có thể thấy ngay S1 = 100 và S2 = 240. Lời giải ban đầu này có tên gọi, “lời giải khả thi cơ bản”, dưới dạng vector có thể diễn đạt: ⎧ X 1 ⎫ ⎧ 0 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪X 2 ⎪ ⎪ 0 ⎪ ⎨ ⎬ = ⎨ ⎬ ⎪ S1 ⎪ ⎪100⎪ ⎪ ⎪ ⎪ ⎪ ⎩ S 2 ⎭ ⎩240⎭ Dưới dạng bảng tính, lợi nhuận được khảo sát theo mẫu sau: Bảng 1.4 Cj 7 5 0 0 X1 X2 S1 S2 Số lượng 0 S1 2 1 1 0 100 0 S2 4 3 0 1 240 - Zj 0 0 0 0 0 (Lợi nhuận) - Cj - Zj 7 5 0 0 - Trong bảng 5.1 các ký hiệu mang ý nghĩa của bài toán đang đặt ra. 22
  24. Cj – Lợi nhuận đơn vị, áp dụng cho cả dòng và cột. Zj – Lợi nhuận. Cj - Zj: Lợi nhuận thuần. Tại bảng 5.1, với giá trị ban đầu như đã nêu, lợi nhuận đang ở mức 0, bài toán chưa đi đến đích. Các động tác tiếp theo phải theo hướng làm cho dòng chứa giá trị lợi nhuận thuần (Cj - Zj ) không âm. Bảng tính chuẩn được lập qua các bước: (1) Chọn biến số tiếp theo để đưa vào cột chuẩn, cột thứ hai tại bảng. Tiêu chuẩn chọn là biến đó có giá trị Cj - Zj lớn nhất. (2) Chọn biến từ cột thứ hai đang ở vị trí ưu tiên phải nhường chỗ cho biến mới tuyển. Để làm điều này, cần tiến hành chia các giá trị trong cột số lượng cho số tương ứng ở cột hai. (3) Tính các giá trị mới cho cột trung tâm vừa thiết lập. (4) Tính các giá trị mới cho các cột còn laiï theo công thức: (Số mới) = (Số cũ) – [ (số tại dòng cũ, trên hoặc dưới số bản lề) x (số tương ứng tại dòng mới, vừa thay trong bước 3) ]. Số bản lề là số nằm tại điểm cắtnhau của dòng trung tâm với cột trung tâm. (5) Tính dòng Zj và Cj - Zj. nếu tất cả giá trị tại dòng Cj - Zj bằng 0 hoặc số âm, bài toán tối ưu đã đến đích. Ngược lại, mọi công việc phải tiến hành lại từ bước đầu tiên. Aùp dụng cách làm vừa nêu vào bảng tính cho phân xưởng gỗ có thể thấy: 1. Vì rằng X1 có giá trị lợi nhuận đơn vị cao nhất, đưa X1 vào cột trung tâm. 2. Chia mỗi số tại cột số lượng với số tương đượng tại cột của X1: 100/2 = 50; 240/4 = 60. Biến X1 thay chỗ S1 tại cột thứ hai. 3. Thay dòng trung tâm bằng cách chia mỗi số trong nó với số chuẩn, nằm tại giao điểm cột trung tâm và dòng trung tâm, cụ thể đây là 2, nơi gặp nhau dòng S1, cột X1: 2/2 =1; ½ =0,5; ½ =0,5; 0/2 =0; 100/2 =50. Dưới dạng bảng các số này sắp xếp như sau: Bảng 1.5 Cj Cột 2 X1 X2 S1 S2 Số lượng 7 X1 1 0,5 0,5 0 50 4. Tính cho dòng S2 mới. 0 = 4 - [4 x 1] 1 = 3 - [4 x 0,5] -2 = 0 - [4 x 0,5] 1 = 1 - [4 x 0] 40 = 240 - [4 x 50] Bảng 1.6 Cj Cột 2 X1 X2 S1 S2 Số lượng 7 X1 1 0,5 0,5 0 50 0 S2 0 1 -2 1 40 23
  25. 5. Tính Zj và Cj - Zj. Bảng 1.7 Cj Cột 2 X1 X2 S1 S2 Số lượng 7 X1 1 0,5 0,5 0 50 0 S2 0 1 -2 1 40 Zj 7 7/2 7/2 0 350 Cj - Zj 0 3/2 -7/2 0 - Tại dòng cuối còn giá trị lớn hơn 0, do vậy lợi nhuận 350x103 chưa phải tối ưu. Cần thiết phải tiếp tục tính lại từ bước thứ nhất cho bảng tiếp theo. 1. Đưa X2 vào cột thứ hai vì 3/2 là số lớn nhất trong dòng cuối. 2. Dòng trung tâm là dòng S2 vì rằng 40/1 = 40 nhỏ hơn 50/(0,5) = 100. 3. Thay thế dòng trung tâm. 4. Tính giá trị mới cho dòng X1. 5. Tính Zj và Cj - Zj. Bảng 1.8 Cj Cột 2 X1 X2 S1 S2 Số lượng 7 X1 1 0 3/2 -1/ 2 30 5 X2 0 1 -2 1 40 Zj 7 5 ½ 3/2 410 Cj - Zj 0 3/2 -1/2 -3/2 - Bài toán qui hoạch trên đây, sau khi giải sẽ đưa ra đáp số X1 = 30 và X2 = 40. Lợi nhuận thu được từ sản xuất tính được bằng VNĐ: 7.103.30 + 5.103.40 = 410.103. Tổng kết quá trình thực hiện trên đây trong khuôn khổ thuật toán của phương pháp simplex. Trước khi tìm max hoặc min của hàm mục tiêu, ví dụ: cTX → maximum cần thực hiện bước chuyển các hạn chế về dạng hệ phương trình tương đương, với sự tham gia của các biến bổ n sung, hay còn gọi là biến lỏng lẻo. Các bất đẳng thức ∑aij x j≤ b i thay bằng các đẳng thức chứa j=1 n thêm si: ∑aij x j+ s i = b i . Bằng cách này có thể xác định bậc của ma trận A tại AX = b bằng m. j=1 Như đã thấy rõ trong ví dụ, ma trận A có thể coi là tập họp hai ma trận dạng [B, N], trong đó B = [ aB1, , aBm ] -ma trận cơ bản và ma trận N kích cỡ m x (n – m). Nếu X là điểm cực trị, bản T T T T thân X có thể viết dưới dạng tương đương như đã thực hiện cho A vừa nêu: X = (XB , XN ) = (b , -1 0), trong đó b = B b ≥ 0. Các biến nằm trong các cột của ma trận B gọi là biến cơ bản, xB1, , xBm như đã biết, các biến khác gọi là biến bổ sung. Từ AX = b; X ≥ 0 có thể chọn giá trị bất kỳ của x, triển khai vecto x dưới dạng như vừa T T T nêu: x = (xB , xN ). Từ AX = b có thể viết BxB + NxN = b. Vecto xB có dạng: -1 -1 xB = B b – B NxN 24
  26. Hàm mục tiêu trở thành: T T T T -1 T T -1 c x = cB xB + cN xN = cB B b + (cN - cB B N)xN = T T T -1 = c X + (cN – cB B N) xN. (1.40) T T -1 Nếu (cN - cB B N) ≥ 0, với xN không âm sẽ thoả mãn bất đẳng thức cTx ≥ cTX và X là điểm cực trị. Phương pháp Simplex áp dụng trong qui hoạch tuyến tính được chương trình hóa. Một trong những minh họa viết bằng ngôn ngữ C được giới thiệu dưới đây, thực hiện theo giải thuật trình bày tại sổ tay tự động hóa tính toán của Wilkinson 3. /* Linear Programming, Simplex method */ #define EPS 1.0e-08 #define FREEALL free_ivector(l3,1,m);free_ivector(l2,1,m);\ free_ivector(l1,1,m); void simplex(a, m, n, m1, m2, m3, icase, izrov, iposv) int m, n, m1, m2, m3, *icase, izrov[], iposv[]; float a; { int i, ip, ir, is, k, kh, kp, m12, nl1, nl2; int *l1, *l2, *l3, *ivector(); float q1, bmax; void simp1(), simp2(), simp(), nerror(), free_ivector(); if ( m != (m1+m2+m3) ) nerror("There is error"); l1=ivector(1, n+1); l2=ivector(1,m); l3=ivector(1,m); nl1=n; for(k=1; k<=n; k++) l1[k]=izrov[k]=k; nl2=m; for (i=1; i<=m; i++) { if (a[i+1][1] < 0.0 ) nerror(" Bad here"); l2[i]= i; iposv[i] = n+1; } for ( i =i; i<=m2; i++) l3[i] = 1; ir=0; if (m2+m3) { ir=1; for (k=1; k<=(n+1); k++) { qi=0.0; for (i=m1+1; i<=m; i++) q1 += a[i+1][k]; a[m+2][k] = -q1; } do { sipm1(a, m+1, l1, nl1, 0, &kp, &bmax); if (bmax <= EPS && a[m+2][1] < -EPS) { *icase = -1; FREEALL return; } else if (bmax <= EPS && a[m+2][1] <= EPS ) { 3 Wilhinson, J.H., Reinsch, C., “Linear Algebra, Vol II of Handbook for Automatic Computation”, Springer-Verlag, 1971. 25
  27. m12=m1+m2+1; if(m12 0.0) goto one; } } } ir=0; m12; if( m1+1 = (n+m1+m2+1) ) { for ( k=1; k = (n+m1+1) ) { kh = iposv[ip] -m1 -n; if (l3[kh] ) { l3[kh] = 0; a[m+2][kp+1] += 1.0; for( i=1; i <= m+2; i++) a[i][kp+1] = -a[i][kp+1]; } } } is=izrov[kp]; izrov[kp] = iposv[ip]; iposv[ip] = is; } while (ir); } for );;) { simp1(a, 0, l1, nl1, 0, &kp, &bmax); if (bmax <= 0.0) { *icase = 0; FREEALL return; } simp2(a, n, l2, nl2, &ip, kp, &q1); if (ip == 0) { *icase =1; FREEALL return; 26
  28. } simp3( a, m, n, ip, kp); is = izrov[kp]; izrov[kp] = iposv[ip]; iposv[ip] = is; } } Những vấn đề thường gặp được giải quyết bằng các phương pháp qui hoạch tuyến tính trong thực tế sản xuất có thể là: - Vận tải hàng hoá, vật tư từ điểm nguồn đến nơi tiêu thụ với chi phí vận chuyển thấp nhất. - Phân công công việc cho người làm thích hợp, giảm chi phí sản xuất đến mức thấp nhất. - Tổ chức sản xuất hợp lý, chọn tuyến đường ngắn nhất trong mạng để thực hiện các công việc trao đổi, phân công. Bài toán vận tải đề cập m nguồn, năng lực sản xuất aj đơn vị ( j = 1, 2, , m); và n đích, năng lực tiêu thụ bk đơn vị nguồn (k = 1, 2, , n). Giá thành vận chuyển đơn vị nguồn từ j đến k là cjk, số lượng vận chuyển xjk từ nguồn j đến đích k. Bài toán vận tải dạng này được đưa về phép lấy minimum của hàm sau: m n Z = ∑∑cjk x jk —>minimum. j==11k trong đó: n ∑ xjk = aj ; j = 1,2, , m k =1 m ∑ xjk = bk ; j = 1,2, , n j=1 và xjk 0 cho tất cả j và k. Trong bài toán này cần thiết phải đặt ra điều kiện: m n ∑aj = ∑ b j j=1 k =1 Bài toán phân công công việc giúp cho việc phân n công việc cho n người hoặc nhóm người thực hiện. Chi phí nhân công của việc thứ j do nhóm k thực hiện là cjk. Bài toán tìm giá nhân công thấp nhất được đưa về dạng: m n Z = ∑∑cjk x jk —>minimum. j==11k trong đó: n ∑ xjk =1; j = 1,2, , m k =1 m ∑ xjk =1; j = 1,2, , n j=1 và xjk 0 hoặc 1, cho tất cả j và k. 27
  29. 1.8 QUI HOẠCH PHI TUYẾN 1.8.1 Hàm một biến Hàm một biến, ví dụ f(x), có minimum cục bộ tại điểm x0 nếu tồn tại đại lượng vô cùng bé δ thoả mãn điều kiện, nếu |x – x0| 0, f”(x ) > 0. (*) Khi đạo hàm bậc hai biểu diễn bằng hai công thức (*) đổi dấu, các điểm vừa chọn sẽ là những điểm mà f(x) mang giá trị maximum. Trường hợp đạo hàm bậc hai vừa nêu bằng 0, bài toán chưa xác định. Để làm rõ hơn vấn đề, có thể xét hàm f(x) khi phân thành chuỗi Taylor. h 2 f()()'() x+ f − fx = hf x + f''( x )+ (1.42) 0 0 0 !2 0 * Nếu tại điểm x0 (hoặc x ) hàm đạt minimum, biểu thức bên trái của (1.42) sẽ âm với giá trị nhỏ bất kỳ của h (|h| 0 trong mọi trường hợp, và nếu f”(x0) > 0, tại x0 f(x) sẽ đạt minimum. Ví dụ: hàm f(x) = x3 – 2x2 + x + 1 đạt cực trị tại các vị trí sau. Tại x = 1/3 đạo hàm bậc hai f’’(x) = 6x – 4 = - 2 0 hàm f(x) đạt minimum. 1.8.2 Hàm nhiều biến Trường hợp biến của hàm không chỉ 1 như phần vừa trình bày mà gồm n biến, ký hiệu x1, x2, x3, , xn, hay dưới dạng vecto sau đây: T x = {x} = [x1, x2, x3, , xn ] hàm nhiều biến có dạng: 28
  30. f(x) = f(x1, x2, x3, , xn ) (1.43) Đạo hàm bậc một của hàm nhiều biến có tên gọi gradient của hàm, là vecto gồm các thành phần đạo hàm riêng ∂f/∂x1, ∂f/∂x 2, ,∂f/∂xn, ký hiệu bằng ∇f(x). Khái niệm đạo hàm bậc hai như đã gặp trong trtường hợp hàm một biến khi áp dụng cho hàm nhiều biến được hiểu là ma trận cỡ (n x n) chứa các thành phần: ∂ 2 f H ij = (1.44) ∂xi ∂ x j được ký hiệu H(x), có tên gọi ma trận Hess. Khai triển thành chuỗi Taylor của hàm nhiều biến đưa đến dạng: n ∂f 1 n n ∂ 2 f f()() x0 + h − f x0 = ∑ hi (x1 , , xn ) + ∑∑hi h j (x1 , , xn )+ i=1 ∂xi !2 i==11 j ∂xi ∂ x j T T = h ∇f(x0) + ½ h H(x0)h + (1.45) Điều kiện cần để hàm nhiều biến có minimum tại x0 là: ∇f(x0) = 0, có nghĩa là ∂f() x 0 =0;(i = 1, , n ) (1.46) ∂xi T trong trường hợp này hiệu giữa f(x0 + h) – f(x0) xác định theo ½h H(x0)h. Nếu H(x0) dương, biểu thức này sẽ dương đối với tất cả h. Do vậy có thể viết điều kiện cần và đủ để hàm nhiều biến kể trên đạt minimum là: ∇f(x0) = 0, H(x0) dương. Điều kiện đạt maximum: ∇f(x0) = 0, H(x0) âm. 2 2 2 Ví dụ: Xác định cực trị hàm ba biến sau f(x) = x1 + x2 + x3 -4x1 -8x2 – 12x3 + 100. Gradient được tìm sau đây: ⎧ 2x1 − 4 ⎫ ⎪ ⎪ ∇f(x) = ⎨ 2x2 − 8 ⎬ = {}0 tại x1 = 2; x2 = 4; x3 = 6. ⎪ ⎪ ⎩2x3 −12⎭ Ma trận Hess: ⎡2 0 0⎤ ⎢ ⎥ H(x) = ⎢0 2 0⎥ là ma trận dương, các giá trị trên đường chéo chính bằng 2. ⎣⎢0 0 2⎦⎥ Do vậy tại điểm (2, 4, 6) hàm f(x) đạt minimum. 29
  31. 1.8.3 Xác định minmum hoặc maximum hàm một biến Phương pháp tìm đoạn vàng trong bài toán một chiều Phương pháp có tên gọi rất hay “Golden Section Search”, trong thực tế là phương pháp chia đoạn có nghiệm (max hoặc min) làm hai phần để tiến dần đến đích, không khác cách làm như đã trình bày về phương pháp bisection, chương hai. Chiều dài đoạn chứa điểm cực trị của hàm được tìm trong: LLLj−1 = j + j+1 (a) Nếu tỷ lệ giữa các phân đoạn cuối là hằng số, thoả mãn điều kiện: LLL j−1 =j = j+1 = = τ (b) LLj j+1 L j+2 chúng ta có thể viết: LL j−1 =1 + j+1 , có nghĩa τ = 1 + 1/τ (c) LLj j Từ đó: τ3 - τ - 1 = 0 và nghiệm của phương trình bậc ba này sẽ là: τ = ( 1+ √5)/2 ≈ 1,61803399. Kết quả cho thấy rằng: L L j−1 =τ2 ; j−2 = τ 3 vv (d) L j+1 L j+1 L Có thể viết như sau 1 = τ n−1 , hay là: Ln L L = 1 (e) n τ n−1 Nếu xác định đúng đoạn có chứa cực trị, ví dụ trong (x0, x3), trong đó xác định được hai giá trị của hàm f(x) tại x1 và x2, với x 1 f2: Phân đoạn được chọn để tiếp tục tìm là (x1, x3 ) x0 x1 x2 x3 Hình 1.2 Chương trình sau minh họa cách tìm điểm cực trị hàm f(x) = -e-xLn(x) theo phương pháp phân đoạn vàng vừa nêu. /* 8/91 for nonlinear programming TCN */ 30
  32. #include #include #include /* phuong trinh chinh*/ float funceval( float x ) { return -exp(-x)* log(x); } void main() { float err, i, inita, initb ; float x0, x1, x2, x3, t1, t2, f1, f2; printf("\n Chuong trinh tinh toi uu- pp. Golden Section "); err = 0.000001; printf("\n Gan 2 gia tri cua han che a, b:="); scanf("%f %f", &inita, &initb); t1= 0.381966; t2= 1.0 - t1; x0 = inita; x1= inita + t1*(initb-inita); x2=inita+t2*(initb-inita); x3=initb; f1= funceval(x1); f2=funceval(x2); printf("\n x1=%8.3f f(x1)= %12.7f x2=%8.3f f(x2)=%12.7f",x1, f1, x2, f2); L1: printf("\n %8.4f %8.4f %8.4f %8.4f %12.7f %12.7f",x0, x3, x1, x2, f1, f2); if (f2 > f1) { i=x2-x0; x3=x2; x2=x1; x1=x0+t1*i; f2=f1; f1= funceval(x1); if ( i > err) goto L1; } else { i=x3-x1; x0=x1; x1=x2; x2=x0+t2*i; f1=f2; f2=funceval(x2); if ( i > err) goto L1; } printf("\n x = %12.7f f(x)= %15.8f", x1, f1); } Nghiệm bài toán được tìm x = 1,762852744 ; f(x) = -9,72601324E-02; Nghiệm chính xác của bài toán là x = 1.76322211. 31
  33. Phương pháp Newton Phương pháp kinh điển tìm vị trí cực trị hàm một biến là xác định dạo hàm bậc nhất và gán cho đạo hàm bằng 0. Giải phương trình vừa xác lập sẽ tìm đúng giá trị của x để f’(x) = 0. Trường hợp chung, giải phương trình dạng này không phải giản đơn. Một trong những cách tìm nghiệm gần đúng của phương trình mới xác lập là phương pháp Newton. PA =tangψ = ϕ '( x ) ( ) a TA 0 Hình 1.3 từ đó: PA ϕ()x TA = = 0 ( ) b ϕ'()x0 ϕ'()x0 ϕ()x0 x1= x 0 − ( ) c ϕ'(x0 ) ϕ()x1 x2= x 1 − ( ) d ϕ'()x1 ϕ()xr và trường hợp tổng quát: xr+1 = x r − ( ) e ϕ'()xr Chương trình tìm giá trị cực tiểu của hàm số f(x) = ½x2 – sinx theo phương pháp Newton được trình bày tiếp theo. Trong chương trình sử dụng các hàm sau: ϕ(x) = f’(x) = x – cosx và ϕ’(x) = 1 + sinx khi tính. 32
  34. /* 8/98 for nonlinear programming */ #include #include #include /* phuong trinh chinh*/ double funceval( double x ) { return 0.5*(x * x) - sin(x); } /* dao ham bac 1 */ double funcd1( double x ) { return x - cos(x); } /* dao ham bac 2 */ double funcd2( double x ) { return 1 + sin(x); } void main() { /* phuong phap Newton, xu ly phi(x) = y'(x) */ double err, init, x, f, d, ff ; printf("\n Chuong trinh tinh toi uu"); err = 0.000001; printf("\n Gan gia tri ban dau:="); scanf("%f", &init); do { x = init; f = funcd1(x); d = funcd2(x); init = x - f/d; printf("%20.8f%20.8f\n", x,init); } while (abs(init - x) > err); x = init; f = funcd1(x); d = funcd2(x); ff = funceval(x); if (d > 0) printf("\n Mimimum tai x = %15.3f", x) ; if (d < 0) printf("\n Maximum tai x = %15.3f", x); printf(" f(x)= %20.8f", ff); } 33
  35. 1.8.4 Phương pháp sử dụng gradient Giả sử vector x của (*) có thể chuyển vị từ x đến x +h.d, trong đó d - hướng chuyển dịch, h - bước chuyển. Thay đổi hàm f(x) sau chuyển vị của x được tính như sau: df = f ( x1 + δx1, x2 + δx2, , xn + δxn) - f ( x1, x1, , xn) ∂f ∂f ∂f = δx1 +δx2 + + δxn (1.47) ∂x1 ∂ x2 ∂xn Trong trường hợp đơn giản nhất phương trình trên được viết thành: df = ∇fxdx() cosθ (5.13) trong đó θ- góc giữa vector ∇f(x) và dx. Trường hợp θ = 180° hướng của dx trùng với -∇f(x). Với bước nhỏ dọc đường đẳng mức, f( x1 + dx1, x2 + dx2, , xn + dxn) = f( x1,x2, ,xn ) biểu thức df có dạng: n ∂f T df = ∑ dj = ∇ f( x ) d = 0 (1.48) j=1 ∂x j Quá trình tìm kiếm cực trị tiến hành như sau: xi+1 = xi - λi∇f(xi), (1.49) trong đó λi - giá trị của λ - hằng số Lagrange, xác định từ quan hệ: ϕ(λ) = f| xi - λi∇f(xi)|, (1.50) Giá trị của λi xác định theo một trong các cách tìm nghiệm hàm một biến đã trình bày. Giải thuật của Powell và Fletcher Bắt đầu Xi,i = 0 di = -g(xi) Tìm λi i= i+1 xi+1 = xi + λidi Min ? x* = x i+1 Hình 1.4 Giải thuật này nhằm tìm hướng đi đển tiến đến mục tiêu nhanh nhất, trong đó sử dụng hàm nguyên thủy và cả gradient của nó lúc khảo sát. Sơ đồ khối cho giải thuật này như tại hình 1.4. 34
  36. Giải thuật Davidon-Fletcher-Powell Trong khuôn khổ phương pháp DFP, các hàm liên tục được viết dưới dạng gần đúng của hàm bậc hai: F(x) = a + xT b + ½ xT H x, với a – hằng số, b – hằng vecto. Hàm đạt cực trị tại x*, khi đó: ∇F(x*) = b + Hx* = 0 (a) Từ đó x* = -H-1 b. (b) Dưới dạng chuỗi Taylor, hàm kiểu này được viết thành: T F(x) = f(x0) + (x – x0) ∇f(x0) + ½ (x – x0) H(x0)(x – x0) (c) Công thức (a) giờ đây có thể viết thành: ∇F(x0) = G(x0)(xm – x0) = 0 (d) -1 Từ đó: xm = x0 - G (x0) ∇F(x0), (e) -1 hoặc dưới dạng: xm = x0 - H (x0) g(x0), với g() thay vào vị trí ∇F(). Nghiệm bài toán tìm dưới dạng: -1 xi+1 = xi - H (xi )g(xi ) (*) -1 hoặc xi+1 = xi - λiG (xi )g(xi ) ( ) Các bước thực hiện như sau: (1) Tại bước thứ j chọn xj và hướng dương của ma trận đối xứng Hj. (2) Chọn dj = -Hjgj làm hướng tìm kiếm cực trị. (3) Thực hiện tìm cực trị theo hướng xj + λjdj và xác định minimum hàm f(xj + λjdj ) để có hàm λj. (4) Đặt: νj = λjdj (5) Đặt xj+1 = xj + νj (6) Tìm f(xj +1) và gj+1. Nếu các đại lượng |g j+1 | hoặc | νj | đủ bé, sẽ dừng công việc. Ngược T lại, tiếp tục thực hiện theo thủ tục trên. Tại đây có thể lưu ý đến quan hệ gj+1 νj = 0. (7) Đặt uj = gj+1 - gj (8) Tính ma trận H tại bước thứ j+1 theo công thức: Hj+1 = Hj + Aj + Bj, trong đó: T T Aj = νjνj / (νj uj); T T Bj = -Hjuj uj Hj / ( uj Hj uj ) (9) Tăng j thêm 1 và quay lại bước 2. Thể hiện phương pháp DFP bằng chương trình bạn đọc có thể tham khảo listing dưới đây. 35
  37. Bắt đầu từ xi, Hi (i = 0) di = Hi gi Xác định λi, minimum hàm F( xi + λi di ) xi+1 = xi + λi di = xi + νi Tìm gi+1 i = i + 1 Từ H,i νi u i=g i+1 –gi để Tạo ma trận Hi+1 Sai Đúng Dừng ⏐ν⏐ #include #include #include #define MAXN 10 double x[MAXN], g[MAXN]; double funceval( double x[] ) { return 0; } double gval( double x[], double g[] ) { return 0; } 36
  38. double DFPoptimization( void ) /* Coded in Hochiminh City 1992 Prepared by Dr.Tran Cong Nghi */ { register int i,j; int cc,n; double H[MAXN][MAXN], p[MAXN], tt[MAXN], d[MAXN], v[MAXN],mm[MAXN]; double y[MAXN],g1, gr, gp, gq, fp, fq, qx, cg, fr,wk,dk, auk, auxz, waux,dd,varh; double w, bb, g2,g3, q[MAXN]; cc=0; printf("\nEnter number of variables N= "); scanf("%d", &n); for( i=0; i 1.0) qx= 1.0; for( i=0; i 1.0) qx =1.0; varh=qx; bb=varh; L7: for(i=0; i<n; i++) x[i]= q[i]= p[i] + bb* d[i]; fq = funceval(x); g2= geval(x,g); cg=0; for( i=0;i<n; i++) cg += g[i]* d[i]; 37
  39. if (gq > 0 || fq >fp ) goto L83; varh = 2* varh; goto L7; L83: auxz = 3*(fp-fq)/varh; auxz += (gp+gq); waux= auxz * auxz - gp*gq; if( waux 00 ) goto L99; varh -= dd; for ( i=0; i < n; i++) p[i] = x[i]; fp= funceval(x); gp= gr; g1= geval(x,g); goto L83; L99: varh = dd; for(i=0; i<n; i++) q[i] = x[i]; L11: auk = dk = wk =0.0; for( i=0; i<n; i++) { tt[i] = g[i] -tt[i]; v[i] = x[i]- y[i]; } for( i=0;i<n; i++) { mm[i] =0.0; for ( j=0; j< n; j++) mm[i] += H[i][j] * tt[j]; auk += mm[i] * tt[i]; wk += v[i] * tt[i]; dk += v[i] * v[i]; } if ( auk == 0.0 || wk == 0.0 ) goto L12; for( i=0; i <n; i++) for(j=0; j <n; j++); H[i][j] += ( -mm[i] * mm[j]/ auk + v[i]*v[j]/ wk); L12: cc++; if( sqrt( dk) < 0.5e-5 || g3 < 0.1e-5) goto L3; goto L4; L3: printf("Iteration number %d\t and f = %20.8f\t",cc,funceval(x) ); for( i=0; i<n; i++) printf("%d%20.8f",i,x[i]); return funceval(x); } 38
  40. 1.8.5. Phương pháp tìm trực tiếp (không qua giai đọan tính gradient) Ý tưởng phương pháp hết sức đơn giản, cố gắng bằng mọi cách thử nghiệm tìm điểm xk, tại đó thỏa mãn điều kiện || xk+1 - xk || < ε. Bước tiến hành đầu tiên của phương pháp là chọn vector d1, d2, , dn, cùng điểm xuất phát x1, và hãy đặt y1 = x1, k = j = 1, sau đó bắt tay vào tính lặp. Gia đoạn tính lặp: (1) Giả sử λj là hằng số Lagrange thỏa mãn lời giải tối ưu f(yj+ λdj). Tìm tiếp giá trị yj+1 = yj + λjdj . Nếu j < n hãy thay j thành j + 1 và quay lại bước đầu. Nếu j = n hãy tiến đến bước (2). (2) Đặt xk+1 = yn+1. Nếu || xk+1 - xk || < ε dừng các phép tính. Trường hợp ngược lại, hãy đặt y1 = xk+1, j=1, thay k thành k+1 và quay về bước (1). Những phương pháp chính trong phần này gồm: - Phương pháp Hooke-Jeeves, - Phương pháp Neldel-Mead, hay còn gọi phương pháp Simplex công bố trong “A Simplex Method for Function Minimisation”, Comp. Jour.,1965. - Phương pháp của Rosenbrock, - Phương pháp của Box hay còn gọi phương pháp Complex. Trong các thủ tục tính nêu trên, khi xử lý những bài toán trong miền hạn chế chúng ta thường gặp khái niệm hằng số Lagrange và khái niệm lồi, lõm hàm đa biến. Trong tài liệu này sẽ không trình bày cách xác định λj cũng như điều kiện Kuhn-Tucker. Về hàm Lagrange và điều kiện của Kuhn- Tucker liên quan miền lồi hàm đa biến được tìm thấy trong các tài liệu chuyên ngành sau: “ A new derivation of the Kuhn-Tucker conditions”, Operations Research, 12, 1964. H.W. Kuhn and A.W. Tucker, “Non linear Programming”, hội thảo khoa học tại trường đại học Berkeley, 1951. Phương pháp Hooke-Jeeves, Phương pháp này được giới thiệu lần đầu năm 1961, còn giá trị đến ngày nay. Quá trình tìm gồm nhiều bước, quanh quẩn điểm cơ bản. Thủ tục tiến hành như sau: (a) Chọn điểm cơ bản b1 và bước hj cho mỗi biến xj, j =1,2,3, , n. (b) Tính f(x) tại điểm cơ bản b1 và xác định hướng tìm kiếm. 1. Tính f(b1), 2. Tính f(b1 +h1e1 ). Nếu động tác này làm cho hàm f() nhỏ hơn, thay b1 bằng b1 +h1e1. Ngược lại, thực hiện phép tính theo f(b1 - h1e1), và nếu giá trị của f() nhỏ hơn sẽ thay b1 thành b1 - h1e1. nếu các động tác trên không mang lại hiệu quả làm nhỏ hơn hàm f(), điểm chuẩn b1 sẽ được giữ nguyên giá trị ban đầu và chuyển hướng sang tìm trên trục x2. Thay vì bước tính f(b1 +h1e1 ) như đã làm ở trên, lần này tính f(b1 + h2e2 ). Cần thiết khảo sát cho tất cả n biến, nếu công việc đòi hỏi, sau đó sẽ chuyển qua điểm cơ bản mới b2. 3. Nếu b2 = b1, có nghĩa là không làm cho f() giảm, phải làm lại từ đầu, cũng quanh quẩn điểm cơ bản ban đầu b1 song bước phải nhỏ hơn. Kinh nghiệm tính cho thấy, nên giảm bước chừng 10 lần. 4. Trường hợp b2 ≠ b1, tiến hành các bước nêu dưới đây quanh điểm cơ bản sau. (c) Sử dụng đầy đủ thông tin thu nhận từ phần trên để tiến hành công việc theo hướng tốt nhất. Thủ tục làm như sau: 39
  41. 1. Theo hướng b2 - b1 để tiếp tục công việc. Tính giá trị hàm f() tại điểm P1 = b1 + 2(b2 - b1). Trường hợp tổng quát tính theo: Trường hợp tổng quát tính theo: Pi = bi + 2(bi+1 - bi). 2. Khảo sát hàm chung quanh điểm P1 (Pi). 3. Thực hiện các động tác nêu trên cho điểm cơ bản b3 nếu cần. (d) Thực hiện các động tác trên cho đến khi bước trở thành quá nhỏ, đạt giới hạn cho phép, nếu cần. Chương trình tính sau đây được soạn để tìm điểm cực tiểu hàm Rosenbroke, theo phương pháp Simplex của Neldr-Mead: /* example for opt. */ #include #include #include #define sqr(X) (X) * (X) #define IMAX 20 int k,m,N; double xb[IMAX], xh[IMAX], xl[IMAX], stp[IMAX]; double fb; double funceval( double xb[] ) { return 100.0*sqr( xb[1] - xb[0]*xb[0] ) + sqr(1.0 - xb[0]) ; } void DirectSearch( int N, double *xb, double fb, double *xh, double *xl, double *stp, int flag, int ifns, int conv ) /* Direct Search procedure in C Coded in Ho Chi Minh City 1992 by Dr. Tran Cong Nghi */ { register int i, j ; int fns, iexp, minfal, rvs[IMAX]; double stor, temp, pat, fbp, fp; double a, b, tol, tol1; double bp[IMAX]; for (i=0; i xh[i] ) xb[i] = xh[i]; if ( xb[i] < xh[i] ) xb[i] = xl[i]; if ( xh[i] == xl[i] ) stp[i] = 0.0; bp[i] = xb[i]; } fns = 1; fb= funceval( xb ); 40
  42. fp = fbp = fb; L2: minfal = 0; if ( fns > ifns ) goto L21; if ( flag > 0 ) { printf ("fns = %d fbp = %20.5lf\n", fns,fbp); for ( i = 0; i xh[i]) || (xb[i] 1 ) { printf("i= %d fb = %20.5lf\n", i, fb); for (j =0; j = fp - tol1 * fabs(fp) ) goto L4; if ( iexp == 0 ) stp[i] *= a; L3: fp = fb; goto L11; L4: xb[i] = stor - stp[i]; if ( (xb[i] > xh[i]) || (xb[i] 1) { printf("i = %d fb = %20.5lf\n", i, fb); for (j=0; j = fp - tol1*fabs(fp) ) goto L5; if (iexp == 0 ) stp[i] *= a; rvs[i] =1; goto L3; L5: xb[i] = stor; if ( iexp == 1) goto L10; stp[i] *= b; temp = fabs(xb[i]/stp[i])* tol; if ( 1.0 = temp ) goto L10; L 8: stp[i] *= temp; goto L10; L 9: temp = 1.0- 16 / fabs( stp[i]); if ( 1.0 == temp ) goto L10; else if ( 1.0 < temp ) goto L8; else goto L11; L10: ++minfal; L11: j = 0; } 41
  43. if ( (fbp-fp) 2) printf("%d", fns); for ( i=0; i xl[i] ) ? xb[i]: xl[i]; xb[i] = ( xb[i] 2 ) { printf("fp = %25.6lf\n", fp); for (i=0; i = N ) goto L18; if (flag > 2 ) printf(" base point exploratory mode failure\n"); for ( i=0; i 2 ) printf("pattern mode exp failure\n"); goto L2; L18: if ( fp 0 ) { printf("The optimum value found\nfp = %25.6lf\n",fp); for (i=0; i 0 ) printf("fns = %d\n",fns); goto L23; L21: conv = 1; if ( flag > 0 ) { printf("number function estimates exceeded %d\n",fns); print f("fbp = %25.6lf\n",fbp); for (i=0; i < N; i++) printf("bp[%d] = %25.6lf\n",i, bp[i]); } fb = fbp; for ( i=0; i < N; i++) xb[i] = bp[i]; L23: return; } 42
  44. main() { int ii,infs,fb,flag,conv; printf("Enter inputs: N, ifns, flag\n"); scanf("%d %d %d", &N, &infs, &flag); printf(" estmates of xb\n"); scanf("%lf %lf %lf", &xb[0], &xb[1], &xb[2]); printf("Enter xhi\n"); scanf("%lf %lf %lf", &xh[0], &xh[1], &xh[2]); printf("Enter xli\n"); scanf("%lf %lf %lf", &xl[0], &xl[1], &xl[2]); for (ii=0; ii < N; ii++) stp[ii] = 0.24 * ( xh[ii] - xl[ii] ); DirectSearch(N, xb, fb, xh, xl, stp, flag, infs, conv); printf("\n\tOutputs:\n"); printf("fb = %20.6lf\n", fb); for (ii=0; ii < N; ii++) printf("xb[%d] = %20.6lf\n", ii, xb[ii]); } 1.8.6. Phương pháp dùng hàm phạt penalty Phương pháp dùng hàm phạt nằm trong phần tính tối ưu (The Sequential Uncostrained Minimisation Technique, viết tắt SUMT) do A.V. Fiacco và G.P. McCornick đề xướng. Để giải hàm f(x) như đã trình bày trên, trong phương pháp SUMT thực hiện phép biến đổi: chuyển z = f(x) thành Z = f(x) + P(x). (1.51) Hàm P(x) được gọi là hàm phạt (Penalty). Hàm P(x) đồng thời mang đặc tính rào chắn, nó có thể còn được gọi là hàm ba-rie (rào chắn). Bây giờ tiến hành tìm cực trị của hàm Z trong miền hạn chế do đầu đề đặt ra, cụ thể hơn trong hàng rào chắn của hàm ba- rie. Hàm P(x) thường đựơc viết dưới dạng: m P(x) = r 1 (1.52) ∑ cx() j=1 j trong đó r - đại lượng mang giá trị dương. Hàm Z = ϕ(x,r) dưới đây có dạng: m Z = ϕ(x,r) = f(x) + r 1 (1.53) ∑ cx() j=1 j Yêu cầu đặt ra cho vector r là, r phải là đại lượng vô cùng nhỏ để ảnh hưởng của hàm P(x) rất nhỏ tại điểm đạt cực trị. Từ đó có thể coi điểm mà hàm ϕ(x,r) không hạn chế đạt cực trị, trùng với điểm cực trị của hàm f(x) cùng các hạn chế. Ví dụ đơn giản nhất về hàm mục tiêu, hàm phạt dưới dạng hàm rào cản được minh họa bằng bài toán đơn thuần, không gắn với hiện tượng vật lý nào, như sau: f(x) = x → min. Các điều kiện hạn chế: x ≥ 2, tương đương với x - 2 ≥ 0. 43
  45. Có thể nhìn thấy ngay lập tức hàm f(x) bậc nhất, khi x tăng, f(x) tăng, do vậy hàm chỉ có thể đạt giá trị nhỏ nhất tại x =2, và f(x) min = 2. Để giải bài toán theo phương pháp hàm phạt, tiến hành xây dựng hàm rào cản P(x). Bài toán được đưa về dạng: m r Z = ϕ(x,r) = f(x) + r 1 = x - ∑ cx() x − 2 j=1 j Lời giải bài toán được minh họa trên hình bên. Trên đồ thị thấy rõ, với r nhỏ dần, từ r= 1, sau đó 0,25, 0,01 điểm Q tiến từ Q1, Q2, và Q3. Giá trị Q đạt được tại x = 2 là điểm gặp nhau của f(x) và đường hạn chế x =2. Phạm vi hạn chế nằm bên phải đường x =2 như thể hiện trên hình. Một trong các cách tìm x tại đó hàm ϕ(x,r) đạt cực trị là khảo cứu x để đạo hàm bậc nhất bằng 0. dϕ r =1 − dx (x − 2) 2 Giá trị để dϕ/dx = 0 chính là x = 2 ± r Hình 1.6 Khi đó: d 2ϕ 2r = dx 2 (x − 2)3 và điểm để hàm đạt minimum là x = 2 + √r, trong miền hạn chế. Giá trị của hàm tại điểm đang xét 2 + 2√2. Ví trí các điểm Qi, i = 1,2,3 xác định trên đồ thị. Tọa độ điểm Q1 sẽ là (3; 4), Q2 tại (2,5 ;3), còn Q3 là (2,1; 2,2). Như vậy khi r → 0, hàm ϕ(x, r) của miền hạn chế tiến đến 2 như đã thấy. Cách tìm điểm cực trị tiến hành theo thứ tự sau. Nếu x1 và x2 là nhữ ng điểm nằm trong vùng được hạn chế, thỏa mãn cj(x1) ≥ 0 và cj(x2) ≥ 0, với i = 1,2, , m, với mọi giá trị của 0 ≤ θ ≤ 1 sẽ thỏa mãn bất đẳng thức sau: cj(θ.x2 + (1 - θ).x1) ≥ θ.cj(x2) + (1 - θ).cj(x1) ≥ 0, cho hàm cj(x) lồi. Và như vậy điểm x2 + (1 - θ).x1 trong phạm vi 0 < θ < 1 cũng sẽ nằm trong phạm vi tìm kiếm. Ngoài ra hàm 1/cj. (x) cũng là hàm lồi cho tất cả giá trị x nếu x thỏa mãn cj(x) ≥ 0. Nếu h(x) = 1/ cj(x) thì 44
  46. −∇cxj () ∇h(x) = 2 (1.54) [()]cxj Ma trận Hess của h(x) có dạng: T Cx() 2∇∇cxcxjj() () H( ) = - + (1.55) x 23 [()]cxj [()]cxj 2 ∂ cxj () trong đó C(x)ik = là hàm Hess của cj(x). ∂∂xxik Nếu ký hiệu p - vector bất kỳ, phương trình chứa p sau đây cũng là phương trình có nghĩa: T T 2 T pCxp() 2[()pcx∇ j ] p H(x)p = + (1.56) 2 3 [()]cxj [()]cxj T trong đó luôn luôn thoả mãn p H(x)p ≥ 0. Với ma trận Hess H(x) dương thì 1/cj(x) phải là lồi trong toàn miền. Giả sử rằng x*1, x*2, , x*n là điểm cực trị của hàm ϕ(x, r) với các giá trị giảm dần r1, r2, ,rk, , cho đến 0. Khi đó chuỗi điểm x*1, x*2, , x*k, , hội tụ về lời giải bài toán tìm cực trị với hạn chế kiểu c đã cho khi rk → 0. Từ đó: lim xk = x* (1.57) và lim[minϕ (x , rk )]= f ( x *) (1.58) rk →0 trong đó x* - điểm hàm f(x) đạt cực trị, với hạn chế đã đặt ra. Trong trường hợp có nghiệm các phép tính đưa về dạng: m 1 f(x* ) → f(x*) và r → 0 (1.59) k k ∑ * j=1 cxjk() 3 Ví dụ:Tìm giá trị nhỏ nhất của hàm f(x1, x2 ) = (1/3) (x1 + 1) + x2, với hạn chế: x1 - 1 ≥ 0; x2 ≥ 0. Tìm hàm ϕ(x, r) dạng (5.19): 1 ⎛ 1 1 ⎞ 3 ⎜ ⎟ ϕ(,)x r=() x1 +1 + x2 + r⎜ + ⎟ 3 ⎝ x1 −1 x2 ⎠ Điều kiện cần để hàm ϕ() đạt minimum sẽ là: 2 r r (x1 + 1) − 2 =0;1 −2 = 0 (x1 + 1) x2 Nghiệm của hai phương trình này: 45
  47. x( r )= 1 + r ; 1 x2 () r= r Giá trị nhỏ nhất của hàm ϕ(x,r) tính từ công thức: 3 3 ⎧1 ⎡ 2/1 ⎤ ⎡ 2/1 ⎤ ⎫ ϕ *()r =⎨ ()1 +r +1 +r + r(1 + r ) +1 ⎬ (1.60) ⎩3 ⎣⎢ ⎦⎥ ⎣⎢ ⎦⎥ ⎭ Khi r → 0 hàm x1(r ) → 1; x2(r) → 0 và ϕ*(r) → f(1; 0) = 8/3. Giải thuật Fiacco và McCornick tóm tắt như sau: m Tìm giá trị nhỏ nhất của hàm ϕ(x,r) = f(x) + r. 1 = f(x) + r. P(x) (5.26) ∑ cx() j=1 j Gradient của hàm ϕ(x,r) tìm dạng sau: ∇ ϕ(x,r) = ∇f(x) + r ∇P(x) (1.61) Từ đó có thể viết: ∇f(x)T∇f(x) + 2r.∇f(x)T∇P(x) + r2∇P(x)T∇P(x) (1.62) và giá trị tối thiểu tìm cho vector r ban đầu: −∇fx()T ∇ Px () r = (1.63) ∇∇Px()T Px () Các giá trị thử nghiệm lần sau có thể là rk+1 = rk / C, với C ví dụ bằng 10. Bắt đầu từ x0 r = r0 Tìm Min. của hàm ϕ(x, r) Tại x * k Ñ Tại xk* là Min ? Dừng k = k +1 1 Sai r = r / c k+1 k x = xk* Hình 1.7 46
  48. Điều có thể nói ở đây để tìm giá trị cực tiểu của hàm ϕ(x,rk+1) tốt nhất nên dùng phương pháp gradient. Cụ thể hơn, trong các phương pháp gradient, với trường hợp vừa nêu phương pháp DFP (Davidon-Fletcher-Powell) thích hợp hơn cả. Chương trình giới thiệu tiếp sau đây nhằm khảo sát hàm: f(x) = (x1 – 1) (x1 – 2) (x1 – 3) + x3 Cùng hạn chế: x1 ≥ 0; x2 ≥ 0; x3 ≥ 0; 2 2 2 x3 - x1 - x2 ≥ 0; 2 2 2 x1 + x2 + x3 - 4 ≥ 0; x3 ≤ 0; ⎛ 1 1 1 1 1 1 ⎞ ϕ(x, r) = f(x) + r⎜ + + + + + ⎟ ⎜ 2 2 2 2 2 2 ⎟ ⎝ x3 − x1 − x2 x1 + x2 + x3 − 4 5 − x3 x 1 x 2 x 3 ⎠ /* The Sequential Unconstrained Minimisation Technique for nonlinear programmming A.V. Fiacco and G.P. McCormick ref. SIAM J. Appl. Math. 17,1969 */ #include #include #include #define M 10 #define N 20 /* double x[N], c[M], cg[N], g[N], r; */ double funceval( double x[], double c[], int m, double r, double a, double b ) { register int i; a = ( x[0] -1) * ( x[0] -2) * ( x[0] -3 ) + x[2]; b =0.0; for ( i=0; i< m; i++) b += 1.0/c[i] ; return ( a+ r*b); } double geval( double x[], double g[], double cg[], double c[], double r,int m) { register int i; double ka, kb, kc, g0; ka = x[0] -1; kb = x[0] - 2; kc = x[0] - 3; g[0] = ka* kb + kb* kc + kc* ka; cg[0] = -(-2*x[0]/c[0]*c[0]) + 2* x[0] /(c[1]*c[1]) + 1.0/ (c[3]*c[3]); g[0] += r* cg[0]; g[1] = 0.0; cg[1] = -(-2*x[1]/c[0]*c[0]) + 2* x[1] /(c[1]*c[1]) + 1.0/ (c[4]*c[4]); /*cg[1] = 0*/ g[1] += r* cg[1]; g[2]= 1.0; cg[2] = -(-2*x[2]/c[0]*c[0]) + 2* x[2] /(c[1]*c[1]) + 1.0/ (c[2]*c[2]); 47
  49. cg[2] -= 1.0/ (c[5] * c[5]); g[2] += r* cg[2]; g0=0; for ( i=0; i 0 ) printf("\nFirst point is not legal. Program terminated"); exit(0); } 48
  50. t = b = r =0.0; cc =0; g1 = geval(x, g, cg, c, r, m ); for( i =0; i 0.0 ) printf( "\nWarning: function is increasing!"); L =2; for(i=0; i = 0 ) goto L73; ic[ii] =1; s++; L /= 1.05; for( i=0; i 0 ) goto L65; hh =L; for(i=0; i <n; i++) x[i] = q[i] = p[i]+ hh*d[i]; for ( ii=0; ii <m; ii++) constrain(x, c, ii ); fq= z = funceval(x, c, m, r, aa, bb ); g2 = g0= geval(x, g, cg, c, r, m); gq =0; for(i=0; i <n; i++) gq += g[i] * d[i]; if( gq < 0.0 && fq < fp ) goto L11; goto L12; L11: 49
  51. for(i=0; i 0.0 ) goto L19; hh -= dd; for( i=0; i < n; i++) p[i] = x[i]; fp =z; gp =gr; g1 =g0; goto L12; L19: hh =dd; for(i=0; i<n; i++) q[i] = x[i]; fq = z; gq = gr; g2 = g0; goto L12; L14: for( i=0; i <n; i++) { u[i] = g[i]- u[i]; v[i] =x[i] -y[i]; } for( i=0; i <n; i++) { mm[i] =0; for( j=0; j < n; j++) mm[i] += h[i][j] * u[j]; kk += mm[i] * u[i]; wk += v[i] * u[i]; } if( kk ==0.0 || wk ==0.0 ) goto L15; for(i=0; i<n; i++) for( j=0; j < n; j++) h[i][j] += -mm[i]* mm[j] /kk + v[i] * v[j] /wk; L15: cc++; if( fabs( (ff - z)/ff < 1.0e-5) )goto L16; ff = z; goto L5; L16: if( r*bb < 1.0e-5 ) goto L18; r /= 10.0; goto L4; L18: for( i=0; i < n; i++) printf("%d %20.8f\n",i, x[i]); printf("\n Value of f(x) = %20.8f",aa); } } 50
  52. Chương 2 TÍNH NỔI VÀ TÍNH ỔN ĐỊNH TÀU 2.1 TÍNH NỔI TÀU THỦY 2.1.1. Kích thước chính và các hệ số thân tàu Phân biệt các tên gọi liên quan đến chiều dài tàu sau: Chiều dài toàn bộ tàu, Lt hoặc Loa, là khoảng cách đo từ điểm xa nhất trước tàu đến điểm xa nhất sau lái. Hình 2.1 Chiều dài đường nước kết cấu LKW, đo trên đường nước thiết kế, kể từ điểm tiếp nước ở mũi tàu đến điểm tiếp nước phía sau lái. Chiều dài giữa hai trụ Lpp, là khoảng cách đo trên mặt đường nước, tính từ trụ lái đến trụ mũi. Trên tàu vỏ thép trụ lái được hiểu là trục đi qua trục quay bánh lái, còn trụ mũi đi qua điểm cắt nhau của đường nước thiết kế với mép ngoài trên lô mũi tàu. Với các tàu có vách đuôi nằm nghiêng so với mặt cơ bản qua đáy (vách T ), trụ lái nhận đi qua đường cắt của vách nghiêng với đường nước thiết kế, tính trên mặt cắt dọc giữa tàu. Chiều rộng tàu lớn nhất Bmax, là khoảng cách lớn nhất đo tại mặt cắt ngang tại khu vực rộng nhất của tàu, tính từ điểm xa nhất bên mạn trái đến điểm xa nhất bên mạn phải của tàu. Hình 2.2 Chiều rộng tàu B, thuật ngữ chuyên ngành bằng tiếng Anh viết đầy đủ là Breadth moulded, là khoảng cách đo từ mạn trái đến mạn phải tàu, tại mặt cắt ngang tàu đi qua mặt rộng nhất của tàu. Với tàu có mặt cắt hình U hoặc V, vị trí đo nằm tại mép boong. Với tàu dạng ω chiều rộng tàu đo tại vị trí rộng nhất của mặt cắt. 51
  53. Chiều cao tàu, ký hiệu bằng D hoặc H, là khoảng cách đo theo chiều thẳng đứng, tính từ mép trong của tấm ki chính đến mép trên của xà ngang boong mạn khô. Với tầu nhiều boong, boong mạn khô được hiểu là boong có kết cấu kín nước, có hệ thống đậy kín các lỗ khoét trên boong và các lỗ khoét bên mạn, nằm ở vị trí cao nhất. Mớn nước tàu ký hiệu bằng d hoặc T, đo trên trục thẳng đứng, tính từ đường cơ bản qua đáy tàu, đến đường nước thiết kế. Với tàu đáy bằng mớn nước tiêu chuẩn đo tại giữa tàu. Phân biệt các tên gọi thường dùng sau. Mớn nước d (chiều chìm), thuật ngữ chuyên ngành trong tiếng Anh gọi là draught moulded (tiếng Mỹ: draft molded) đo từ đường cơ bản. Chiều cao đo từ mép dưới sống chính gọi là keel draft, còn mớn nước trung bình dm là giá trị trung bình cộng của mớn nước đo tại trụ lái và mớn nước đo tại trụ mũi. Mớn nước lái đo tại trụ lái, tính cả chiều nghiêng của sống chính, nếu có. Mớn nước mũi đo tại trụ mũi, tính cả độ nghiêng của sống chính. Mạn khô Chiều cao mạn khô tàu là hiệu số giữa chiều cao và mớn nước tàu: Fb = D - d hoặc H - T Hệ số đầy Quan hệ giữa kích thước chính của tàu với thể tích phần chìm, diện tích đường nước, diện tích mặt giữa tàu vv đựợc thể hiện qua các hệ số đầy. Hệ số đầy đường nước Hình 2.3 Hình 2.4 CW hoặc α, là tỉ lệ giữa diện tích mặt đường nước được vỏ tàu giới hạn và diện tích hình chữ nhật có cạnh là chiều dài và chiều rộng đường nước. Nếu ký hiệu AW - diện tích mặt đường nước, L - chiều dài tàu, đo tại đường nước, B - chiều rọng tàu, hệ số CW tính theo công thức: V C = B LBT×× Hình 2.5 52
  54. AW CW = LxB Hệ số đầy sườn giữa tàu, CM hoặc β, là tỉ lệ giữa diện tích phần chìm của sườn giữa tàu AM với diện tích hình chữ nhật ngoại tiếp nó, cạnh BxT AM CM = BxT Hệ số đầy thể tích, CB hoặc δ, là tỉ lệ giữa thể tích phần chìm của tàu V với thể tích hình hộp ngoại tiếp nó. Hệ số CB tính theo công thức: V CB = LxBxT Hệ số đầy lăng trụ, CP hoặc ϕ, là tỉ lệ giữa thể tích phần chìm tàu V so với ống trụ dài bằng chiều dài đường nước L, diện tích mặt trụ AM. V CB CP = hay là CP = AM xL CM Hệ số đầy trụ đứng, CV hoặc χ, là tỉ lệ giữa thể tích phần chìm so với trụ đứng cao T, mặt trụ AW. V CB CV = hay là CP = AxT W CW Đường hình vỏ tàu Đường hình lý thuyết của vỏ tàu được biểu diễn trong hệ toạ độ gắn liền với vỏ tàu như trên hình. Trục OZ hướng lên trên. Trục OX trùng với chiều dọc tàu, hướng về trước, còn trục OY hướng sang mạn phải. Tâm của hệ toạ độ đặt tại giao điểm ba mặt phẳng: mặt cắt ngang giữa tàu, mặt cắt dọc giữa tàu và mặt cơ bản qua đáy tàu. Hình 2.6 Đường lý thuyết miêu tả vỏ tàu được qui ước vẽ trong bản vẽ hai chiều 2D, bao gồm các phần sau. Hình chiếu các vết cắt dọc tàu do mặt phẳng dọc giữa tàu và các mặt phẳng song song với mặt này tạo thành. Cụm vết cắt này nằm phía trái, trên. Hình chiếu các vết cắt vỏ tàu qua các đường nước, nằm phía trái, dưới. Hình chiếu các mặt cắt ngang tàu, gọi là các sườn lý thuyết, nằm phía phải, trên. Tính các đại lượng hình học vỏ tàu 53
  55. Từ đường hình lý thuyết tiến hành tính các giá trị đặc trưng hình học vỏ tàu. Thứ tự tính toán chia làm hai giai đoạn: (1) các đại lượng đặc trưng trong mặt đường nước và (2) các đại lượng đặc trưng trong sườn tàu. Sau hai phần tính vừa nêu tiến hành tính toán cho toàn tàu. Đại lượng hình học đường nước Biểu diễn đường nước bất kỳ của tàu dưới dạng đường cong dạng y = f(x), các phép tính đại lượng hình học đường nước được đưa về dạng sau. Hình 2.7 Diện tích đường nước AW. b AW = ydx ∫a Tấn trên 1cm chiều chìm (Tons per cm Immersion – TPcm) A Trong môi trường nước ngọt: TPcm = W 100 1025, × A Trong môi trường nước biển: TPcm = W 100 2 trong đó: AW tính bằng m . b Moment tĩnh so với trục Oy: moy = xydx ∫a b 1 2 Moment tĩnh so với trục Ox: mox = y dx ∫a 2 Tọa độ tâm diện tích đường nước, tính đến trục Oy: b ∫ xydx LCF = a b ydx ∫a Tọa độ tâm diện tích đường nước, tính đến trục Ox: 1 b ∫ y2dx TCF = 2 a b ydx ∫a Momen quán tính so với trục dọc: b 2 IxL = 2 ydx ∫a 54
  56. Momen quán tính mặt đường nước so với trục O’y’ cách Oy một đoạn LCF tính theo công thức trên sẽ là: 2 IL’ = IL -(LCF) .AW Momen quán tính mặt đường nước so với trục ngang: 2 b 3 IyT = dx 3 ∫a Trong các biểu thức trên y mang giá trị ½ chiều rộng vỏ tàu tại vị trí đang xét. 2.1.2 Tỉ lệ Bonjean Diện tích mặt sườn tính đến mớn nước Z. z Sz()= 2 ydz . ∫0 Momen tĩnh so với trục Oy của mặt sườn: z mz()= 2 yzdz . ∫0 Tâm diện tích mặt sườn thuộc phần chìm đến mớn nước Z tính theo công thức: z mz() ∫ yzdz. VC() z ==0 z Sz() y.dz ∫0 Với mỗi sườn tàu, từ kết quả tính diện tích phần chìm và momen tĩnh phần chìm so với đáy, có thể vẽ hai đường cong miêu tả biến thiên của hai giá trị trên theo chiều chìm Z. Tập hợp toàn bộ các đường cong kiểu này, lập cho tất cả sườn tính toán sẽ được đồ thị có tên gọi tỉ lệ Bonjean. 2.1.3 Thể tích phần chìm và các đại lượng liên quan đến thể tích Tính thể tích phần chìm được tiến hành theo một trong hai cách: 1- Tính từ dưới lên trên cơ sở dữ liệu của tất cả đường nước hoặc 2- Tính theo chiều dọc tàu, sử dụng dữ liệu các sườn làm cơ sở. Hình 2.9 Thể tích phần chìm, tính đến mớn nước Z: 55
  57. z V()zAz= W ().dz ∫0 trong đó V- thể tích phần chìm, AW(z) - diện tích đường nước. Sử dụng tỉ lệ Bonjean khi tính thể tích phần chìm: L/ 2 V()zSx= ().dx ∫−L/ 2 trong đó S(x) diện tích sườn, tính đến mớn nước Z, đọc từ biểu đồ Bonjean. Momen thể tích phần chìm so với mặt phẳng qua đáy tàu: z MAzXOY= W () zdz ∫0 Momen thể tích phần chìm so với mặt giữa tàu: L/ 2 MSxYOZ = () xdx ∫−L/ 2 Tọa độ tâm nổi phần chìm: Chiều cao tâm nổi, ký hiệu KB hoặc ZB: z ∫ AzzdzW () Mz() KB≡= Z 0 =XOY B z Vz() AzdzW (). ∫0 Hoành độ tâm nổi, ký hiệu LCB hoặc XB: L ∫ Sx() xdx Mz() LCB≡= X 0 =YOZ B L Sx(). dx Vz() ∫0 Các đường cong tính nổi Kết quả tính các đặc trưng hình học vỏ tàu được tập họp trong một bảng vẽ chung mang tên gọi các đường cong tính nổi của tàu. Thuật ngữ chuyên ngành để chỉ đồ thị dạng này không giống nhau ở nhiều nước, trong đó có nước ta. Một nước có nền công nghiệp đóng tàu phát triển, có ảnh hưởng lớn đến sách vở của Việt nam gọi đây là các đường cong yếu tố đường hình. Trong tài liệu chính thức của tổ chứchàng hải quốc tế IMO, họ đường cong này mang tên gọi bằng tiếng Anh hydrostatic curves, có nghĩa các đường thủy tĩnh của tàu. Các đường cong bắt buộc có mặt trong đồ thị này: 1. Thể tích phần chìm ∇ hoặc V. 2. Lượng chiếm nước Δ = γ∇ hoặc ký hiệu cách khác D = γ.V 3. Hoành độ tâm nổi LCB. I 4. Chiều cao tâm nghiêng KM tính theo công thức: KM=+ KB BM =+ KB T V 5. Hoành độ tâm đường nước LCF IL 6. Momen làm nghiêng dọc tàu 1 cm, ký hiệu MTRIM hoặc MTcm: MTcm =γ L 56
  58. 7. Chiều cao tâm nổi KB I 8. Chiều cao tâm nghiêng dọc: KM=+ KB BM =+ KB L LLV 9. Các hệ số đầy (béo) thân tàu. Hình 2.10 Các đường thủy tĩnh Hình 2.11Đồ thị Bonjean 57
  59. 2.1.4 TÍNH CÁC ĐƯỜNG THỦY TĨNH TRÊN MÁY CÁ NHÂN Các đại lượng hình học trình bày trên được chia làm ba nhóm khác nhau lúc tính. (1) Tính diện tích, momen tĩnh, momen quán tính, hệ số đầy đường nước, momen chúi đơn vị vv trong mỗi đường nước. (2) Tính diện tích phần chìm các sườn, momen tĩnh so với đáy, so với mặt giữa tàu cho mỗi sườn, thực hiện trong mặt sườn (3) Tính thể tích phần chìm và các đại lượng liên quan đến thể tích. Các phép tích phân được phân vào hai dạng, tích phân giới hạn xác định dọc chiều dài tầu và tích phân giới hạn trên thay đổi thùy thuộc mớn nước tính toán. Chuẩn bị dữ liệu. 1/ Chọn số sườn tính toán, số đường nước cần thiết khi tính, 2/ Vị trí các sườn tính toán ghi trong hệ toạ độ tương đối, đơn vị tính dL =Lpp/(10 hoặc 20), ví dụ: Bảng 2.4 Thứ tự 1 2 3 NS-1 NS Vị trí sườn #0 #½ #1 #1½ #2 #19½ 20 3/ Vị trí các đường nước theo đơn vị tính dT, ví dụ: Bảng 2.5 Thứ tự 1 2 3 NW Vị trí đường nước 0 ½ 1 4/ Toạ độ vòm đuôi so với trụ lái, ghi lại dưới dạng bảng Bảng 2.6 Thứ tự 1 2 3 NW+1 Toạ độ vòm đuôi 5/ Chiều dài thật của tất cả đường nước tính toán, 6/ Chiều cao của tất cả các sườn tính toán, 7/ Toạ độ vỏ tàu, xác định tại tất cả các sườn tính toán, qua tất cả đường nước tính toán. Toạ độvỏ tàu ( giá trị ½ chiều rộng tàu) ghi dưới dạng ma trận như ví dụ sau: Bảng 2.8 Thứ tự đường nước 0 1/2 1 NW Boong Sườn # 0 #½ # sườn cuối Thứ tự ghi dữ liệu như minh hoạ trên hình trên đây. Chương trình tính thực hiện các phép tính theo thứ tự sau: 58
  60. Tích phân trong mặt đường nước thứ j, j =1,2, , NW Diện tích: AW = 2 ydx (a) ∫L Tâm đường nước: x y dx LCF = ∫L (b) ydx ∫L Momen quán tính dọc, qua trục trung hòa: 2 2 IL = 2 xydx. - a .AW (c) ∫L Momen quán tính ngang: 2 It = ydx3 . (d) 3 ∫L Momen chúi tàu 1 m. MTRIM = IL / L (e) Tích phân trong mặt sườn thứ i, i = 1, 2, , NS Diện tích phần chìm: Z S(z) = 2 ∫ yd. z (f) 0 Momen tĩnh so với đáy: Z MB(z) = 2 ∫ yzdz (g1) 0 Momen tĩnh so với mặt giữa tàu: Z MO(z) = 2 ∫ xydz (g2) 0 Tích phân theo thể tích phần chìm từ 0 đến Z: Thể tích phần chìm: Z V()() z= ∫ Aw z dz (h) 0 Chiều cao tâm nổi: Z ∫ Aw( z ). z . dz KB() z = 0 (i) V() z Hoành độ tâm nổi: 59
  61. Z ∫ Aw( z ). a ( z ). dz LCB() z = 0 (k) V() z Bán kính tâm nghiêng ngang: BM = It / V (l) Bán kính tâm nghiêng dọc: BML = IL /V (m) Các hệ số đầy: CW = AW(z) / L.B (n) CM = SO (z) /T.B (p) CB = V/ L.B.T (q) CP = CB/ CM (r) CV = CB/CW Như đã chuẩn bị tại phần đầu chương này, các phương pháp tích phân số được sử dụng để thực hiện các phép tính trong sơ đồ. Theo cách làm vừa nêu, các phép tích phân thực hiện trong mỗi đường nước thuộc dạng sau: B I = ∫ f() x dx , A Trong đó A vị trí đuôi tàu thuộc đường nước, đo tại mặt đối xứng dọc giữa tàu, B là vị trí lô mũi tàu, đo tại đường nước tính toán, trong mặt đối xứng dọc. Vị trí của A và B thay đổi từ đường nước tới đường nước trên các tàu thông dụng. Các phép tính thực hiện tại mỗi sườn tàu thuộc dạng: z I(z) = ∫ f() z dz , trong đó z biến miêu tả chiều chìm, tính đến đường nước. Tích phân dạng 0 này cần thực hiện theo chương trình. Một trong những hàm viết bằng ngôn ngữ C, giải tích phân dạng vừa nêu được trình bày tại phần đầu chương. 2.1.5 Biểu đồ mang tên Firsov Tên gọi này chỉ thịnh hành tại đất nước đã sinh ra nhà khoa học tàu thủy đáng nể này. Trong sách báo của nước ta tên gọi biểu đồ Firsov được công nhận một cách chính thức, giống như được gọi ở Nga vậy. Tên gọi này không mấy khi xuất hiện trên sách chuyên môn các nước ngoài Nga. Biểu đồ Firsov giúp cho người đọc tìm được thể tích phần chìm hoặc lượng chiếm nước, toạ độ tâm nổi phần chìm cho các trạng thái nghiêng dọc tàu. Nói theo cách dễ hiểu hơn, khi đọc được chiều chìm tàu tại mũi và lái của tàu, người đọc sử dụng biểu đồ Firsov để tìm giá trị thực của lượng chiếm nước D, cao độ tâm nổi KB và hoành độ tâm nổi XB của tàu trong trạng thái ấy. Biểu đồ Firsov trình bày quan hệ giữa thể tích phần chìm hoặc lượng chiếm nước, tọa độ tâm nổi phần chìm cho các trạng thái nghiêng dọc tàu. Khi đọc được chiều chìm tàu tại mũi và lái của tàu, từ biểu đồ Firsov tìm giá trị thực của lượng chiếm nước D, cao độ tâm nổi KB và hoành độ tâm nổi LCB của tàu trong trạng thái ấy. 60
  62. Thủ tục lập biểu đồ Firsov theo thứ tự sau: 1) Xác định mớn nước lái Tl và mớn nước mũi Tm 2) Xác lập đường nước qua hai vị trí trên 3) Sử dụng biểu đồ Bonjean tính thể tích phần chìm và tọa độ tâm nổi theo các công thức: Hình2.13 Biểu đồ Firsov tàu vận tải Thể tích phần chìm của tàu: L V = ∫ axdx() 0 Momen thể tích phần chìm so với đáy: L Mmxd= ()x b ∫ 0 Momen thể tích phần chìm so với mặt cắt ngang giữa tàu: +L/ 2 Mxax⊗ = ∫ .()dx −L/ 2 Chiều cao tâm nổi, so với mặt đáy: L mxdx() M ∫ KB ==b 0 V L ∫ axdx() 0 Hoành độ tâm nổi, tính từ mặt cắt ngang giữa tàu: 61
  63. +L/ 2 xa.() x dx M ∫ LCB ==⊗ −L/ 2 V L ∫ axdx() 0 2.2 ỔN ĐỊNH TÀU Công việc liên quan tính toán ổn định tàu bao gồm: (a) Tính các đại lượng đặc trưng ổn định ban đầu. (b) Tính họ đường cross curves (pantokaren). (c) Xác lập đồ thị ổn định tĩnh, ổn định động. (d) Kiểm tra tính ổn định tàu theo tiêu chuẩn ổn định. 2.2.1 Ổn định ban đầu Hình2.14 Ổn định ban đầu Tàu nghiêng ở góc nhỏ, khoảng cách BB' xác định như sau: v × h BB' = ∇ trong đó v – thể tích nêm tam giác; ∇ - thể tích phần chìm thân tàu. Momen phục hồi: MGph =Δ× Z trong đó GZ - đoạn thẳng nối từ trọng tâm tàu G đến đường tác động lực nổi tại thời điểm xem xét, qua điểm B’. GZ= GM sin ϕ GM là chiều cao tâm nghiêng, M – tâm nghiêng. GZ gọi là tay đòn ổn định. IT Bán kính nghiêng ngang: BMT = V 62
  64. IL ' Bán kính nghiêng dọc: BM L = V Chiều cao của điểm M so với mặt chuẩn qua đáy: KM=+ KB BM hoặc KM = KB + BM Chiều cao tâm nghiêng ngang ban đầu: GM=− KM KG hoặc GM = KB + BM - KG trong đó KG - chiều cao trọng tâm so với mặt chuẩn qua đáy tàu. Tâm nổi B di chuyển trên cung gần cung tròn, bán kính r = BM , tâm tại M. GM=+ KB BM −= KG BM −() KG − KB =− r a trong đó a = KG− KB . GZ = (r - a) sinϕ = rsinϕ - asinϕ Hình 2.15 Tay đòn ổn định tàu GZ và Hình 2.16 chiều cao tâm nghiêng GM 2.2.2 Ổn định tại góc nghiêng lớn Từ đồ thị miêu tả quĩ đạo tâm ổn định M và tâm nổi B trong quá trình tàu nghiêng có thể thấy rõ, tàu nghiêng đến góc đủ lớn, khoảng từ 10 -15° trở lên, tâm M không còn nằm trên trục đối xứng, còn B di chuyển không phải trên cung gần tròn như ban đầu mà theo đường cong không thành luật. Độ tăng tay đòn momen ngẫu lực giữa lực nổi và trọng lực không còn tuyến tính với góc nghiêng mà sang hẵn giai đoạn phi tuyến. Tại các góc nghiêng lớn bán kính ổn định theo nghĩa là khoảng cách theo chiều đứng giữa B’M’, trong đó B’, M’ là vị trí nhất thời ứng với góc nghiêng đang xét, tính theo tỉ lệ giữa momen quán tính đường nước nghiêng và thể tích phần chìm tàu. BM = Iϕ / ∇ Toạ độ tâm nổi B dời đến vị trí mới: ϕ Y = ∫ BM()cosϕϕ dϕ 0 ϕ Z = ∫ BM()sinϕϕ dϕ 0 63
  65. Khoảng cách GMϕ khi tàu nghiêng tính theo công thức đã trình bày tại phần trên: d GMϕ = GZ = - Ysinϕ + (Z - KB) cosϕ- a.cosϕ dϕ Momen ổn định tĩnh: Δ×GZ =Δ×( BR − BT ) Từ quan hệ: v ×=∇×hh' BR có thể tính: v × hh' BR = ; BT= BG sin ϕ ∇ Momen ổn định tĩnh tính theo công thức: ⎛⎞v × hh' MB=Δ⎜⎟ −Gsin ϕ ⎝⎠∇ Hình 2.17 Xác định tay đòn ổn định tại góc nghiêng lớn Công thức tính tay đòn ổn định: GZ=−=− BR BT BR BG sin ϕ Họ đường BR=∇ϕ f (,) gọi chung là cross curves hoặc pantokaren. Đồ thị tay đòn ổn định tĩnh GZ = f(ϕ) có dạng: Hình 2.18 Tay đòn ổn định tĩnh GZ = f(ϕ) Momen phục hồi là tích số của lượng chiếm nước D với tay đòn GZ: Mph = D. GZ 64
  66. 2.2.3 Đồ thị ổn định Đồ thị ổn định được hiểu là đường momen ngẫu lực phụ thuộc vào góc nghiêng tàu. Momen này là tích số của tay đòn GZ, biến thiên theo góc nghiêng tàu và lượng chiếm nước, có giá trị không đổi cho một trạng thái khai thác. Do vậy, đồ thị tay đòn GZ = f(ϕ) cũng được coi là đồ thị ổn định. Trong các phần tiếp theo của tài liệu khái niệm đồ thị ổn định được dùng cho cả hai đồ thị vừa nêu. Momen ngẫu lực vừa nêu được gọi là momen phục hồi. Mỗi đồ thị ổn định có thể thấy rõ ba đoạn. Đoạn đầu, phụ thuộc vào chiều cao tâm ổn định ban đầu, GM với giá trị dương, đồ thị sẽ tăng trưởng nhanh hoặc chạâm tùy thuộc vào độ lớn GM . Miền tăng trưởng kết thúc tại vị trí cực trị GZmax ứng với góc ϕm. Trong đoạn tiếp theo đường cong vẫn mang giá trị dương song theo xu hướng giảm dần. Giai đoạn này kết thúc tại góc lạên ϕV, tại đây GZ trở lại giá trị 0. Đoạn tiếp theo GZ < 0. Trên đồ thị ổn định có thể biểu diễn độ dốc đường ổn định dưới dạng tiếp tuyến với đồ thị ổn định, xuất phát của tiếp tuyến từ gốc tọa độ. Chiều cao tâm ổn định ban đầu đọc trên đường tiếp tuyến tại vị trí ϕ = 57,3 Giá trị lớn nhất tay đòn ổn định GZmax . Góc ứng với tay đòn lớn nhất ϕm . Góc lặn ϕV Hình 2.19 Tay đòn ổn định tĩnh 2.3 THUẬT TOÁN XÁC LẬP HỌ ĐƯỜNG CROSS CURVES (PANTOKAREN) Để xác định tọa độ tâm nổi B’ cho góc nghiêng bất kỳ, xét trong hệ tọa độ chung toàn tàu, tiến hành cách rời rạc hóa bài toán theo các bước sau. 1/ Phân chia toàn bộ chiều dài tàu thành những phân đoạn, có chiều dài phân đoạn ngắn hơn nhiều lần chiều dài tàu. Chiều dài các phân đoạn không nhất thiết bằng nhau. Mỗi phân đoạn rất ngắn kiểu này được coi như một khối trụ dài đúng bằng chiều dài phân đoạn, mặt cắt ngang của lăng trụ đúng như mặt cắt ngang giữa lăng trụ. 2/ Trong mỗi phân đoạn tiến hành tính thể tích phần chìm, tâm nổi phần chìm so với đáy, so với mặt ngang chuẩn, cụ thể so với mặt cắt ngang giữa tàu. Thuật toán tính thể tích và momen tĩnh lăng trụ. 65
  67. Hình 2.20 Tại mỗi mặt sườn tiến hành kẻ nhiều đường nước nghiêng dưới góc ϕ so với mặt đáy, cắt sườn tàu. Tại mỗi chiều chìm Z, tính trên trục OZ, kẻ đường nước song song với mặt thoáng nước tĩnh. Xác định các giá trị bổ trợ a,b,c theo công thức: z a =−t (a) tgϕ 1 by=+(a) (b) 2 cy=−a (c) trong đó: y - nửa chiều rộng tàu, đo tại mớn nước z, tại sườn đang xét; t - khoảng cách từ giao điểm mặt đường nước nghiêng tại đáy đến mặt cắt dọc giữa tàu. Diện tích mặt sườn, phần nằm dưới đường nước nghiêng: T A()xcz= ()dz (d) ∫0 Momen tĩnh so với mặt đáy: T MxB()= czzd () z (e) ∫0 Momen tĩnh so với mặt cắt ngang giữa tàu: T MxC ()= czbzd ().()z (f) ∫0 3- Tính thể tích phần chìm tàu và tọa độ tâm nổi của phần chìm tàu, nằm dưới một đường nước nghiêng, khi đã xác định A(x), MB(x), Mc(x). Thể tích phần chìm: ∫ MxdxC () ∫ MxdxB() V = Axd()x; RB = L ; SB = L (g) ∫ A()xdx A()xdx L ∫ ∫ L L Tay đòn hình dáng Lk , hình 2.20 tính theo công thức: Lk = KN sinφ=SBsin φ + RBcos φ (h) 66
  68. Mọi giá trị c(z) đo trên bản vẽ phải là những giá trị thật, có nghĩa c(z) ≥ 0. Những trường hợp thường gặp khi đọc c(z) và cách hiệu chỉnh như sau: Nếu c ≥ 2y giá trị thật của c = 2y; b = 0 Nếu c ≤ 0 thì: c = 0 và b = 0; Hình 2.21 Những trường hợp thường gặp khi xác định c Sơ đồ tính được giới thiệu tại hình 2.22. Hình 2.22 Xây dựng họ đường cross curves tàu nhiều thân Họ đường cross curves LK = f(V, ϕ) lập cho trường hợp Vi; i = 1,2, với góc nghiêng ϕ (hoặc còn ký hiệu ) thay đổi từ 0 đến góc bất kỳ, ví dụ đến 90°. 67
  69. Tay đòn hình dáng LK được đo từ điểm K giao điểm của mặt cắt dọc giữa tàu với đường cơ bản nằm ngang, đến hướng tác động lực qua tâm nổi B’ của phần chìm tàu trong thời điểm tính toán, ứng với góc nghiêng cho trước (H.2.23). Điểm K cố định trong mọi trường hợp tính toán. Hình 2.23 Tại mỗi mặt sườn tiến hành kẻ nhiều đường nước nghiêng dưới góc Φ so với mặt đáy, cắt sườn tàu. Tại mỗi chiều chìm Z, tính trên trục OZ, kẻ đường nước song song với mặt thoáng nước tĩnh. Xác định các a, b, c theo công thức trình bày tại (a), (b), (c) cho tàu một thân. Nhóm công thức này áp dụng cho “tàu ảo” nằm phía trong của tàu hai thân, giúp tính “tay đòn ảo”. Các đại lượng nhận được sau khi tính được ký hiệu a’, b’, c’. Diện tích mặt sườn, phần nằm dưới đường nước nghiêng: T A(x) = czdz() ∫0 Momen tĩnh so với mặt đáy: T MB(x) = cz() zdz ∫0 Momen tĩnh so với mặt cắt ngang giữa tàu: T Mc(x) = cz().() bzdz ∫0 Tính thể tích phần chìm tàu và tọa độ tâm nổi của phần chìm tàu, nằm dưới một đường nước nghiêng, khi đã xác định A(x), MB(x), Mc(x). Tay đòn hình dáng LK tính theo công thức: Lk = KN sinφ=SBsin φ + RBcos φ C R O S S C U R V E S Phi 10.00 20.00 30.00 volume lever volume lever volume lever (m3) (m) (m3) (m) (m3) (m) 225.207 0.492 238.849 0.807 243.231 1.070 68
  70. 196.814 0.541 219.026 0.896 220.268 1.193 165.752 0.561 195.195 0.989 192.002 1.342 135.570 0.589 168.254 1.078 159.623 1.503 106.524 0.627 139.006 1.155 124.880 1.655 79.444 0.687 109.489 1.222 89.778 1.760 54.515 0.775 82.233 1.306 57.490 1.833 33.650 0.916 58.244 1.386 30.868 1.958 16.712 1.016 37.892 1.477 12.127 2.211 3.845 1.163 21.191 1.618 2.291 2.475 VN 272.135 0.322 256.417 0.716 224.447 1.171 Phi 70.00 80.00 90.00 volume lever volume lever volume lever (m3) (m) (m3) (m) (m3) (m) 238.902 1.756 235.783 1.795 232.471 1.769 216.962 1.815 212.260 1.824 208.762 1.763 193.416 1.880 187.070 1.859 183.360 1.764 168.553 1.950 160.436 1.902 156.440 1.775 142.483 2.029 132.384 1.961 128.462 1.802 115.358 2.123 103.897 2.036 101.062 1.835 88.216 2.234 76.790 2.112 75.185 1.867 62.260 2.335 51.201 2.180 50.825 1.903 37.789 2.395 27.441 2.205 28.083 1.936 16.913 2.361 7.630 2.132 8.362 2.022 VN 100.745 2.182 79.235 2.105 60.994 1.889 Đồ thị pantokaren (Cross Curves) có dạng. Hình 2.24 Cross Curves (Pantokaren) Dựng đồ thị ổn định trên cơ sở pantokaren Đồ thị ổn định được dưng dưới dạng đường GZ = f(ϕ); Momen ổn định dựng dưới dạng momen phục hồi bằng tích số của GZ(ϕ) với lượng chiềm nước D = const. 69
  71. GZ = Lk - KG sinϕ Ưùng với mỗi trường hợp V = D/γ = const, từ đồ thị pantokaren dễ dàng đo được Lk, tính theo góc nghiêng, ví dụ 10°, 20°, 30°, Thay giá trị Lk vừa đo được vào biểu thức cuối cùng chúng ta xác lập được dãy giá trị GZϕ, tính cho ϕ = 10°, 20°, 30°, Hình 2.25 Xây dưng đồ thị ổn định tĩnh 2.4 GIỚI THIỆU CHƯƠNG TRÌNH TÍNH TÍNH NỔI TÀU Qui ước đánh số sườn, đường nước. + Sườn #0 qua trụ lái, sườn thứ 10 ( hoặc sườn 20) qua trụ mũi + WL 0 bắt đầu từ đường cơ bản, các WL tiếp theo lên cao dần. + Khoảng cách vòm lái, lấy giá trị + nếu nằm trước trụ lái Chuẩn bị dữ liệu: + NS- tổng số sườn tính toán + NW- tổng số đường nước tính toán. Mặt boong tàu được coi là đường số NW+1 + Số thứ tự các đường nước, số thực, ví dụ:0, 0.5, 0.75, 1, 1.5,2, vv Hình 2.32 + Số thứ tự đường nước, số thực, ví dụ: 0, 0.5, 1, 1.5, 2 vv + Chiều dài thật các đường nước vừa đăng ký + Tọa độ vòm đuôi,tính cho các đường nước vừa đăng ký, dấu theo H.4. + Chiều cao H, cho tất cả sườn tính toán. + Tọa độ vỏ tàu, lấy 1/2 chiều rông tàu. Thứ tự đọc: bắt đầu từ sườn đầu tiên, ví dụ #0; trong mỗi sườn đọc giá trị y của đường nước đầu tiên, ví dụ WL0, đến đường nước cuối cùng là NW+1 ( mép boong). Sử dụng chương trình Nhập dữ liệu ten tau: tên gọi của tàu 70
  72. Tong so suon: số sườn tính toán So duong nuoc: số đường nước tính toán Chieu dai Lpp: chiêu dài giữa hai trụ, m Khoảng suon dL: khoảng sườn lý thuyết, m Khoang duong nuoc dT: khoảng đường nước,m Nhap so thu tu ( số sườn) phai tinh: nhập đủ số chỉ vị trí sườn. Các giá trị cách nhau khoảng trống. Nhap chieu cao ( số sườn): toàn bộ chiều cao của các sườn. Các giá trị cách nhau khoảng trống. Nhap so thu tu (số đường nước): số thứ tự của WL, từ dưới lên boong Chieu dai that ( số đường nước): giá trị chiều dài thật các đường nước, m. Các giá trị cách nhau khoảng trống. Khoang cach vom lai den truc sau ( số đường nước): Khoảng cách cùng giá trị (+), (-) như hình 4 trên. Toa do vo tau: Tọa độ vỏ tàu ( ½ chiều rộng),nhập theo hướng dẫn ghi tại màn hình đối thoại. Toa do diem vao nuoc: Toạ độ diểm để nước tràn vào tàu, giá trị y và z, đo bằng m Ghi chú: Điểm chọn có thể như hình bên cạnh hoặc tại mép boong. Các file chứa dữ liệu: H.DAT - dữ liệu đầu vào, nhập từ bàn phím. Gọi và sửa bằng trình EDIT của DOS. H.OUT - dữ liệu đầu ra STAB.OUT - Bảng tính bằng số tay đòn ổn định tĩnh và ổn định động, tính bằng , cách nhau 10° một lần tính. Giá trị trong bảng tương ứng đồ thị được hiện trên màn hình. Các file sau đây chứa giá trị trung gian FIRSOV.OUT - biểu đồ Firsov CROSS.dat - giá trị pantokaren BON.dat giá trị biểu đồ Bonjean dùng khi vẽ đồ thị. Giá trị tính toán và ký hiệu trên các đồ thị có nghĩa như sau: KB- chiều cao tâm nổi BM- bán kính ổn tâm GM - chiều cao tâm ổn định Tay đòn trong đồ thị pantokaren được xác định như trên hình bên. Lk là khoảng cách từ K đến đường vuông góc với đường nước nghiêng, đi qua tâm nổi phần chìm - điểm B’. Chương trình viết bằng ngôn ngữ C. Dưới đây sẽ trình bày nội dung chương trình tính bằng ngôn ngữ C. 1/ IntegrationV() - tính tích phân có giơí hạn trên thay đổi I(z) = ∫ f(z)dz 2/ Integ() - tích phân giới hạn dưới từ a đến giới hạn trên b: I = ∫f(x)dx 3/ Lag() - Nội suy (ngoại suy) Lagrange 71
  73. 4/ Inputfiles - ghi dữ liệu vào file H.dat 5/ Crossfiles - ghi dữ liệu trung gian vào file cross.dat 6/ FirsoffDiagram - ghi giá trị tính biểu đồ Firsov ra file firsoff.out 7/ BonjeanGr - biểu đồ Bonjean 8/ CrossCurveGr - đồ thị pantokaren 9/ StabilityCurvesGr - đồ thi ổn định tĩnh và động 10/ Stability - Ổn định tàu. Trích đoạn chương trình bằng ngôn ngữ C có dạng như sau. /* Program HydrostaticCurves; */ #include #include #include #include #include #include #include #include #define sqr(A) (A) * (A) char ship[20]; int Nwplus1, Nw, Ns,Index, Midship ; int MaxX, MaxY, MaxColors; float Lpp, dt, dl, yfa, zfa, Cend, Afj, Haftl, Breadth, BxL, d; float Phi, SinPhi, CosPhi, TangPhi, db, dz, an, bn, cn; float yy[9][2]; float aux[100]; float y, sa, sm; float *WL, *Wa, *dis, *VOL, *xb, *KB, *BML, *CW, *CB, *CM, *CP; float *ILc, *It, *BMt, *Af, *Low, *St, *H, *xl, *xa, *zz; float *bb, *cc, *ff, *z, *Aw, *amm, *am, *displa, *trim; FILE *inf, *outf, *cross, *bon ; void StabilityCurvesGr() { char buffer[20]; int k, j; int Nn,x0,y0,x1,y1,x2,y2,y01,y10,ErrorCode,graphdriver,GraphMode ; int p0,KK ; float vmax,lmax,c4, p,p3,p2,p03 ; graphdriver = DETECT; initgraph(&graphdriver, &GraphMode, " \\BGI"); ErrorCode = graphresult(); if ( ErrorCode != 0 ) { fprintf(stderr,"%s",grapherrormsg(ErrorCode)); exit(1); } MaxY = getmaxy(); MaxX = getmaxx(); MaxColors = getmaxcolor() + 1; vmax =0; lmax =0; for (j =1 ; j lmax) lmax = cc[j]; 72
  74. if (CB[j] > vmax ) vmax = CB[j]; } if (lmax y0) y10=y0; if ( y2 > y0 ) y2=y0; line(j*x1,y10,(j+1)*x1,y2); y10= y01- (int) (p3*p03); y2= y01 - (int) (p3* p2); if ( y10 y0) y10=y0; if ( y2 > y0 ) y2=y0; line(j*x1,y10,x1*(j+1),y2); 73
  75. delay(100); } getch(); closegraph(); } void Stability() { int n,km, i, j, k ; int KK; printf("\nNumber of States N = "); scanf("%d", &n); CB = calloc( 10, sizeof( float )); Wa = calloc( n, sizeof( float )); dis = calloc( n, sizeof( float )); cc = calloc( 10, sizeof( float )); xb = calloc(10, sizeof( float )); xa = calloc(10, sizeof( float )); sa = dim2(9, 11); sm = dim2(9, 11); printf("\nNhap %d gia tri cua Vol: ", n); for (i =0 ; i sa[k-1][0] ) xb[k] =sm[k-1][0]; else { KK=-1; do { KK++; } while ( bn <= sa[k-1][KK] ); db=(-sa[k-1][KK-1] +bn); dz =(- sm[k-1][KK]+sm[k-1][KK-1]); xb[k] = db*dz/( sa[k-1][KK-1]- sa[k-1][KK])+sm[k-1][KK-1]; } fprintf(stdout,"\n\tTRUONG HOP SO %d",i); fprintf(stdout,"\n\tV = %10.2f m3 Zg = %6.2f m",bn,cn ); fprintf(stdout,"\n deg " ); for (k=2 ; k <= 10; k++) fprintf(stdout,"%7d", (k-1)*10); fprintf(stdout,"\n Lk " ); for (k =1 ; k < 10; k++) fprintf(stdout,"%7.2f", xb[k] ); fprintf(stdout,"\n Lst " ); cc[0]=0.0; for (k=1 ; k < 10; k++) { Phi=k*0.174533; SinPhi= sin(Phi); cc[k]= xb[k] - cn*SinPhi; fprintf(stdout,"%7.2f", cc[k]); 74
  76. } integrationv(10,xa,cc,CB); fprintf(stdout,"\n Ld " ); for (k=1 ; k hoac ==>:"); ch = getche(); ch = toupper(ch); if ( ch == '1' ) { printf("\n\tSo suon phai tinh Ns = "); scanf("%d", &Ns); printf("\n\tSo duong nuoc WL = "); scanf("%d", &Nw); location(); printf("\n\tLpp, dL, dT ==>: "); scanf("%f%f%f", &Lpp, &dl, &dt); printf("\n\n\tNHAP SO TT %d SUON PHAI TINH\n", Ns); for (k=0; k < Ns ;k++) scanf("%f", &St[k]); printf("\n\tNHAP CHIEU CAO %d SUON\n", Ns); for (k=0 ; k < Ns; k++) scanf("%f", &H[k]); printf("\n\tNHAP SO TT %d DUONG NUOC\n", Nw); for (j=0; j < Nw; j++) scanf("%f", &WL[j]); printf("\n\tCHIEU DAI THAT %d DUONG NUOC\n", Nw); for (j=0; j < Nw ; j++) scanf("%f", &Low[j]); printf("\n\tKHOANG CACH VOM LAI DEN TRU SAU %d DG NUOC\n", Nw); for (j=0; j < Nw; j++) scanf("%f", &Af[j]); textcolor( 3 ); clrscr(); gotoxy(10,1); printf("\n\t\t\tTOA DO VO TAU\n"); printf("\n\tBAT DAU DOC TU SUON %5.2f DEN SUON %5.2f ",St[0],St[Ns-1]); printf("\nTai moi suon NHAP %d gia tri 0,5B tu day len\n",Nwplus1); for (k=0 ; k < Ns; k++) { printf("\nfor Section no %d = Sation St[%5.2f]\n", k, St[k]); 75