叠片圆形铁芯最大面积问题(趣味问题)

lin2009

铁芯直径为100mm, 由三级硅钢片叠制而成, 求最大截面时的各级片宽和叠厚尺寸(如有多种方案应全部保留)。要求铁芯片宽为5mm的倍数, 末级片宽 >= 30mm。

lin2009

数学模型：
设从里层到外层的铁芯叠片的片宽和厚度分别为(2*x1, 2*y1)1片放置在最中间, (2*x2, y2)2片两侧对称放置, (2*x3, y3)2片对称放置末级。
则由题设条件可得 : 面积为 S = 4*(x1*y1 + x2*y2 + x3*y3) ;
约束条件为 x1 > x2 > x3 >= 15;
            x1, x2, x3 为5的倍数;
            因为是实际物理尺寸, 所以 y1, y2, y3 > 0 且
            sqrt(x1^2 + y1^2) = 50;
            sqrt(x2^2 + (y1+y2)^2) = 50;
            sqrt(x3^2 + (y1 + y2 + y3)^2) = 50; (圆形铁芯直径为100mm)
求 S 的最大值。

这是个简单的优化问题，直接给出参考答案：最大面积为6680.268887 mm^2

20wangz

2# lin2009
Title "Let me try 一 try!";
MultiRun = 50;
Parameter  k(3)[3,10,0],y(3)[0,];
ConstStr  x1=5*k1,x2=5*k2,x3=5*k3;
MaxFunction  4*(x1*y1+x2*y2+x3*y3);
             sqrt(x1^2+y1^2)=50;
             sqrt(x2^2+(y1+y2)^2)=50;
             sqrt(x3^2+(y1+y2+y3)^2)=50;
             k1>k2>k3;


"let me try 一 try!"

====== 结果 ======

重复计算最佳目标函数: 6680.26887381314
重复计算求得最佳目标函数次数: 38
重复计算最差目标函数: 6130.12701937958
重复计算平均目标函数: 6643.37137795577
重复计算总用时: 0:00:19
优化算法: 通用全局优化法(UGO1)
函数表达式: 4*(((5*k1))*y1+((5*k2))*y2+((5*k3))*y3)
目标函数值(最大): 6680.26887439063
k1: 9
k2: 7
k3: 4
y1: 21.7944947282139
y2: 13.9126474203854
y3: 10.1186148057272

约束函数
   1: sqrt(((5*k1))^2+y1^2)-(50) = 4.581444557E-9
   2: sqrt(((5*k2))^2+(y1+y2)^2)-(50) = 4.202753701E-9
   3: sqrt(((5*k3))^2+(y1+y2+y3)^2)-(50) = 4.370029671E-9
   4: k2-(k1) = -2
   5: k2-(k3) = 3

====== 计算结束 ======

lin2009

3# 20wangz

这也是我最初的做法，不谋而合，呵呵。
可以看出都要运行很多次才能得出更大值。


1stopt程序和结果(算法为自动选择)（方式1）
MultiRun = 20;
IntParameter k1 = [ 3, 10 ], k2 = [ 3, 10 ], k3 = [ 3, 10 ];
ConstStr x1 = 5*k1, x2 = 5*k2, x3 = 5*k3;
ParameterDomain = [ 0, 100 ];
MaxFunction 4*(x1*y1 + x2*y2 + x3*y3) ;
k1 > k2;
k2 > k3;
sqrt(x1^2 + y1^2) = 50;
sqrt(x2^2 + (y1+y2)^2) = 50;
sqrt(x3^2 + (y1 + y2 + y3)^2) = 50;
= == == = 结果 = == == =
重复计算最佳目标函数 : 6680.26887322856
重复计算求得最佳目标函数次数 : 3
重复计算最差目标函数 : 6000.04000000567
重复计算平均目标函数 : 6428.89190315853
重复计算总用时 : 0:00:24
优化算法 : 鲁棒简面体爬山法 + 通用全局优化法(SM3)
函数表达式 : 4*((5*k1)*y1 + (5*k2)*y2 + (5*k3)*y3)
目标函数值(最大) : 6680.26887323566
k1:9
k2:7
k3:4
y1:21.7944947177034
y2:13.9126474250109
y3:10.1186148068442
约束函数
   1:k1 - (k2) = 2
   2:k2 - (k3) = 3
   3:sqrt((5*k1)^2 + y1^2) - (50) = 0
   4:sqrt((5*k2)^2 + (y1+y2)^2) - (50) = 0
   5:sqrt((5*k3)^2 + (y1 + y2 + y3)^2) - (50) = 7.105427358E-15
= == == = 计算结束 = == == =

lin2009

或者：
1stopt程序和结果(算法为自动选择)（方式2）

MultiRun = 20;
//  IntParameter x1, x2, x3, y1, y2, y3;
//  Constant x1 = 30, x2 = 40, x3 = 25;
ParameterDomain = [ 0, 100 ];
MaxFunction 4*(x1*y1 + x2*y2 + x3*y3) ;
x1 > x2;
x2 > x3;
x3 >= 15;
round(x1/5) = x1/5;
round(x2/5) = x2/5;
round(x3/5) = x3/5;
sqrt(x1^2 + y1^2) = 50;
sqrt(x2^2 + (y1+y2)^2) = 50;
sqrt(x3^2 + (y1 + y2 + y3)^2) = 50;

//    ====== 结果 ======
//  
//  重复计算最佳目标函数 : 6630.1924933389
//  重复计算求得最佳目标函数次数 : 1
//  重复计算最差目标函数 : 5406.41045463005
//  重复计算平均目标函数 : 5968.0462734463
//  重复计算总用时 : 0:00:28
//  优化算法 : 最大继承法(MIO1)
//  函数表达式 : 4*(x1*y1 + x2*y2 + x3*y3)
//  目标函数值(最大) : 6630.1924933389
//  x1:45
//  y1:21.7944947177034
//  x2:35
//  y2:13.9126474250109
//  x3:25
//  y3:7.59412804650769
//  
//  约束函数
//  1:x1 - (x2) = 10
//  2:x2 - (x3) = 10
//  3:x3 - (15) = 10
//  4:round(x1/5) - (x1/5) = 0
//  5:round(x2/5) - (x2/5) = 0
//  6:round(x3/5) - (x3/5) = 0
//  7:sqrt(x1^2 + y1^2) - (50) = 0
//  8:sqrt(x2^2 + (y1+y2)^2) - (50) = 0
//  9:sqrt(x3^2 + (y1 + y2 + y3)^2) - (50) = 0
//  
//   = == == = 计算结束 = == == =

第一次没有找出最大值！

//    ====== 结果 ======
//  
//  重复计算最佳目标函数 : 6680.26887323433
//  重复计算求得最佳目标函数次数 : 1
//  重复计算最差目标函数 : 5235.93989433958
//  重复计算平均目标函数 : 6125.85245915718
//  重复计算总用时 : 0:01:34
//  优化算法 : 最大继承法(MIO1)
//  函数表达式 : 4*(x1*y1 + x2*y2 + x3*y3)
//  目标函数值(最大) : 6680.26887323566
//  x1:45
//  y1:21.7944947177034
//  x2:35
//  y2:13.9126474250109
//  x3:20
//  y3:10.1186148068442
//  
//  约束函数
//  1:x1 - (x2) = 10
//  2:x2 - (x3) = 15
//  3:x3 - (15) = 5
//  4:round(x1/5) - (x1/5) = 0
//  5:round(x2/5) - (x2/5) = 0
//  6:round(x3/5) - (x3/5) = 1.33226763E-15
//  7:sqrt(x1^2 + y1^2) - (50) = 0
//  8:sqrt(x2^2 + (y1+y2)^2) - (50) = 0
//  9:sqrt(x3^2 + (y1 + y2 + y3)^2) - (50) = 0
//  
//   = == == = 计算结束 = == == =

第二次找到的值比第一次大, 即更优解, 是否是最大最大值呢, 如何判断呢？

lin2009

可能将y1,y2,y3的表达式直接代入优化目标函数（而不是作为约束条件）进行优化会更高效。
实际上这类问题用穷举法更高效，更能确保找到的是最大值。
1＃的参考答案就是用穷举法得出的。

shamohu

用下面代码，可以极短时间内100%得到最好结果：

1. 
   
2. Algorithm = UGO1[100];
   
3. IntParameter k1 = [ 3, 10 ], k2 = [ 3, 10 ], k3 = [ 3, 10 ];
   
4. ConstStr x1 = 5*k1, x2 = 5*k2, x3 = 5*k3;
   
5. ParameterDomain = [ 0, 100 ];
   
6. MaxFunction 4*(x1*y1 + x2*y2 + x3*y3) ;
   
7. k1 > k2;
   
8. k2 > k3;
   
9. sqrt(x1^2 + y1^2) = 50;
   
10. sqrt(x2^2 + (y1+y2)^2) = 50;
    
11. sqrt(x3^2 + (y1 + y2 + y3)^2) = 50;
    

复制代码

迭代数: 29
计算用时(时:分:秒:毫秒): 00:00:00:313
计算中止原因: 达到收敛判定标准
优化算法: 通用全局优化法(UGO1)
函数表达式: 4*(((5*k1))*y1 + ((5*k2))*y2 + ((5*k3))*y3)
目标函数值(最大): 6680.26887438761
k1: 9
k2: 7
k3: 4
y1: 21.7944947280725
y2: 13.9126474206086
y3: 10.1186148056169

约束函数:
   1: k1-(k2) = 2
   2: k2-(k3) = 3
   3: sqrt(((5*k1))^2 + y1^2)-(50) = 4.519804975E-9
   4: sqrt(((5*k2))^2 + (y1+y2)^2)-(50) = 4.261188735E-9
   5: sqrt(((5*k3))^2 + (y1 + y2 + y3)^2)-(50) = 4.343966964E-9

lin2009

用穷举法解这样的题目应该是最简单。总共只有7*6 *5 = 210种情况。
将y1,y2,y3的约束条件代入目标函数，用穷举法的结果和程序如下：

相应的matlab程序如下
运行结果 :
Ok
Elapsed time is 0.000172 seconds.
ans =

           45           35           20       21.794       13.913       10.119       6680.3
% B =
%
%            25           20           15       43.301       2.5245       1.8712       4644.4
%            30           20           15           40       5.8258       1.8712       5378.3
%            35           20           15       35.707       10.119       1.8712       5920.8
%            40           20           15           30       15.826       1.8712       6178.3
%            40           25           15           30       13.301       4.3957       6393.9
%            40           30           15           30           10        7.697       6461.8
%            45           30           15       21.794       18.206        7.697       6569.5
%            45           35           15       21.794       13.913        11.99       6590.2
%            45           35           20       21.794       13.913       10.119       6680.3
           
相应程序 :

1. 
   
2. clear all;
   
3. clc;
   
4. 
   
5. tic;
   
6. smax = 0;
   
7. k = 1;
   
8. ind = zeros(1,7);
   
9. for x3 = 15 : 5 : 45
   
10.     for x2 = (x3+5) : 5 : 45
    
11.         for x1 = (x2+5) : 5 : 45
    
12.             y1 = sqrt(50^2 - x1^2);
    
13.             y2 = sqrt(50^2 - x2^2) - y1;
    
14.             y3 = sqrt(50^2 - x3^2) - y1 - y2;
    
15.             S  = 4*(x1*y1 + x2*y2 + x3*y3);
    
16.             if S >= smax
    
17.                 ind(k,:) = [ x1, x2, x3, y1, y2, y3, S ];
    
18.                 smax = S;
    
19.                 k = k + 1;
    
20.             end
    
21.         end
    
22.     end
    
23. end
    
24. disp('Ok')
    
25. toc
    
26. B = sortrows(ind, 7);  %排序
    
27. format short g;
    
28. B(end,:)                    % 输出最大值
    

复制代码

用穷举的方法得到的确定是最大值，并且可以肯定得出只有一种方案使得面积最大！

20wangz

7# shamohu
运行您的代码，设置Multi=50，能达到最大值的次数只在45次左右，很难达到100%。即使将并行数改到200也不行。

在1stOpt中如何选择算法是个让人挠头的问题，不知道shamohu版主有什么心得，还请不吝赐教。

mrlimail

学习了^_^....