單純形法求解最函式極值問題 matlab程式碼
阿新 • • 發佈:2019-01-13
最近整理以前的程式碼,將以前老師上課的作業程式碼重新整理,分享出來,作業的內容是編寫單純形法,對測試函式進行尋優(極大值或者極小值)。
首先介紹一下單純形法:將上課的ppt轉化為圖片。ppt藍色背景,眼睛快看瞎了
按照ppt的描述編寫演算法如下:
clear all; clc; % mode可以選擇測試函式 % mode = 'exp_test'; % mode = 'Schaffer'; mode = 'Rosen_Brock'; if strcmp(mode,'Schaffer') x=-4:0.1:4; y=-4:0.1:4; [X,Y]=meshgrid(x,y); Z = 0.5-((sin(sqrt(X.^2+Y.^2)).^2)-0.5)./(1+0.001.*(X.^2+Y.^2)).^2; figure(1); surf(X,Y,Z); title('Schaffer Function'); xlabel('X-軸'); ylabel('Y-軸'); zlabel('Z-軸'); end if strcmp(mode,'exp_test') x=-4:0.1:4; y=-4:0.1:4; [X,Y]=meshgrid(x,y); % Z = 1400-(20*(X.^2-Y.^2).^2-(1-Y).^2-3*(1+Y).^2+0.3); % Z = 5*sin(X.*Y)+X.^2+Y.^2; Z = X.*exp(-X.^2-Y.^2)+1; figure(1); surf(X,Y,Z); title('exp_test Function'); xlabel('X-軸'); ylabel('Y-軸'); zlabel('Z-軸'); end if strcmp(mode,'Rosen_Brock') x=-2.048:0.1:2.048; y=-2.048:0.1:2.048; [X,Y]=meshgrid(x,y); % Z = 1400-(20*(X.^2-Y.^2).^2-(1-Y).^2-3*(1+Y).^2+0.3); % Z = 5*sin(X.*Y)+X.^2+Y.^2; % Z = X.*exp(-X.^2-Y.^2)+1; Z = 100.*(Y-X.^2).^2+(1-X).^2; figure(1); surf(X,Y,Z); title('Rosen Brock valley Function'); xlabel('X-軸'); ylabel('Y-軸'); zlabel('Z-軸'); end if strcmp(mode,'Rosen_Brock') variable_min_x = -2.048; % 變數的上下限 variable_max_x = 2.048; variable_min_y = -2.048; % 變數的上下限 variable_max_y = 2.048; else variable_min_x = -4; % 變數的上下限 variable_max_x = 4; variable_min_y = -4; % 變數的上下限 variable_max_y = 4; end % variable_min_x = -4; % 變數的上下限 % variable_max_x = 4; % variable_min_y = -4; % 變數的上下限 % variable_max_y = 4; yita = 0.5; alpha = 1.3; precision = 10^(-3); % Schaffer function % x=variable_min_x:0.1:variable_max_x; % y=variable_min_y:0.1:variable_max_y; % [X,Y]=meshgrid(x,y); % Z = 0.5-((sin(sqrt(X.^2+Y.^2)).^2)-0.5)./(1+0.001.*(X.^2+Y.^2)).^2; % surf(X,Y,Z); % 複合型法尋找schafer functiong 最大值 %% 產生初始的頂點 sige % M=4; % N=3; x=zeros(4,2); for i=1:4 x_1 = variable_min_x+rand*(variable_max_x-variable_min_x); x_2 = variable_min_y+rand*(variable_max_y-variable_min_y); x(i,:)=[x_1,x_2]; end max_iteraton = 250; % 最大的迭代次數 %% % 判斷索取點是否滿足要求: index_error = 1; error = []; optimization_path = []; for gen=1:max_iteraton % 判斷點是否在可行域內 alpha = 1.3; for i=1:4 if ((x(i,1)<variable_min_x)||(x(i,1)>variable_max_x))||((x(i,2)<variable_min_y)||(x(i,2)>variable_max_y)) error(index_error) = i; % 不滿足的點 index_error = index_error+1; end end % error中存放不滿足條件的點的索引值 if isempty(error) == 0 % error 非空 % 對不滿足的點進行處理 norm = []; % 滿足約束條件的點 abnorm = []; % 不滿足約束條件的點 %%%%%%%%%%%----- len_error = length(error); for i=1:len_error abnorm(i,:) = x(error(i),:); % 壞點 end kk = 1; for m = 1:4 %kk = kk+1; for n = 1:len_error if x(m,:) == abnorm(n,:) break; end if n >= len_error norm(kk,:) = x(m,:); kk = kk+1; end end end % 將不滿足條件的點進行調整 % 正常點的形星 norm_center = 0; for m = 1:size(norm,1) norm_center = norm_center + norm(m,:); end norm_center = norm_center/(size(norm,1)); for j = 1:size(abnorm,1) abnorm_new(j,:) = norm_center + yita*(abnorm(j,:)-norm_center); end % 組成新的點 x = [norm;abnorm_new]; end % 計算函式值 for k = 1:4 f(k) = test_function(x(k,:),mode); end [f_min,index_1] = min(f); % 尋找最好的點 [f_max,index_2] = max(f); % 最差的點X_h 次差的點X_g 最好的點X_l X_h = x(index_1,:); X_l = x(index_2,:); cnt = 1; left_points = []; left_two_points = []; for i=1:4 if (i == index_1)||(i == index_2) continue; else left_points(cnt) = f(i); % 除過最好的點和最差的點 left_two_points(cnt,:) = x(i,:); cnt = cnt + 1; end end [left_points_min,index_3] = min(left_points); X_g = left_two_points(index_3,:); % 求形星: cnt = 1; for i=1:4 if i==index_1 continue; else points_except_worst(cnt,:) = x(i,:); cnt = cnt + 1; end end % 計算形星 sum = 0; for i = 1:length(points_except_worst) sum = sum + points_except_worst(i,:); end X_c = sum/length(points_except_worst); % 判斷形星是否在可行域內 if ((X_c(1)>variable_min_x)&&(X_c(1)<variable_max_x))&&((X_c(2)>variable_min_y)&&(X_c(2)<variable_max_y)) % 在可行域內 % 求反射點 X_r = X_c + alpha*(X_c - X_h); while 1 if ((X_r(1)>variable_min_x)&&(X_r(1)<variable_max_x))&&((X_r(2)>variable_min_y)&&(X_r(2)<variable_max_y)) % 反射點在可行域內 break; else if alpha < precision break; end alpha = alpha*0.5; X_r = X_c + alpha*(X_c - X_h); end end else % 形星不在可行域內 % 重新確定上下界 if X_l(1) <= X_c(1) variable_min_x = X_l(1); variable_max_x = X_c(1); else variable_min_x = X_c(1); variable_max_x = X_1(1); end if X_l(2) <= X_c(2) variable_min_y = X_l(2); variable_max_y = X_c(2); else variable_min_y = X_c(2); variable_max_y = X_l(2); end % 產生新的複合型 for i=1:4 x_1 = variable_min_x+rand*(variable_max_x-variable_min_x); x_2 = variable_min_y+rand*(variable_max_y-variable_min_y); x(i,:)=[x_1,x_2]; end for k = 1:4 f(k) = Schaffer(x(k,:)); end [f_min,index_1] = min(f); % 最差的點 [f_max,index_2] = max(f); % 尋找最好的點 % 最差的點X_h 次差的點X_g 最好的點X_l X_h = x(index_1,:); X_l = x(index_2,:); cnt = 1; left_two_points = []; for i=1:4 if (i == index_1)||(i == index_2) continue; else left_points(cnt) = f(i); % 除過最好的點和最差的點 left_two_points(cnt,:) = x(i,:); cnt = cnt + 1; end end [left_points_min,index_3] = min(left_points); X_g = left_two_points(index_3,:); % 求形星: cnt = 1; for i=1:4 if i==index_1 continue; else points_except_worst(cnt,:) = x(i,:); cnt = cnt + 1; end end % 由新的點產生的形星 sum = 0; for i = 1:length(points_except_worst) sum = sum + points_except_worst(i,:); end X_c = sum/length(points_except_worst); % 求反射點 X_r = X_c + alpha*(X_c - X_h); while 1 if ((X_r(1)>variable_min_x)&&(X_r(1)<variable_max_x))&&((X_r(2)>variable_min_y)&&(X_r(2)<variable_max_y)) % 反射點在可行域內 break; else if alpha < precision break; end alpha = alpha*0.5; X_r = X_c + alpha*(X_c - X_h); end end % 反射點計算結束 end % x X-r X_h X_l X_g f_X_r = test_function(X_r,mode); flag_1 = 0; if f_X_r > f_min x(index_1,:) = X_r; % 替換最差的點 else % 重新計算反射點 while test_function(X_r,mode) < f_min alpha = alpha*0.5; X_r = X_c + alpha*(X_c - X_h); if alpha < precision flag_1 = 1; break; end end % 重新替換 if flag_1 == 0 % 新的反射點滿足約束條件 x(index_1,:) = X_r; else % alpha 減小到足夠的範圍任然沒有找到合適的點 % 使複合型向最好的點收縮 再重新計算一次新的複合型 %Nop(); %disp('nono'); x(index_1,:) = X_g; end end % x 更新完畢 如果想畫成凸多邊形,要處理x figure(2); for i=1:4-1 plot([x(i,1),x(i+1,1)], [x(i,2),x(i+1,2)]); hold on; end plot([x(4,1),x(1,1)],[x(4,2),x(1,2)]); hold off; pause(0.01); disp([gen, alpha]); function_value(gen) = test_function(X_l,mode); max_value = test_function(X_l,mode); optimization_path(gen,:) = X_l; end figure(3); plot(function_value); title(['函式值變化曲線',' 最大值:',num2str(max_value)]); xlabel('迭代次數'); ylabel('函式值'); figure(4); plot(optimization_path(:,1),optimization_path(:,2)); title('尋優軌跡'); xlabel('最優點-X'); ylabel('最優點-Y'); % function f = Schaffer(buf) % f = 0.5-((sin(sqrt(buf(1).^2+buf(2).^2)).^2)-0.5)./(1+0.001.*(buf(1).^2+buf(2).^2)).^2; % end % function f = test_function(buf,mode) % if strcmp(mode,'Schaffer') % f = 0.5-((sin(sqrt(buf(1).^2+buf(2).^2)).^2)-0.5)./(1+0.001.*(buf(1).^2+buf(2).^2)).^2; % end % % % if mode=='Rastrigin' % % f = (buf(1).^2-10*cos(2*pi*buf(1))+10) + (buf(2).^2-10*cos(2.*pi.*buf(2))+10); % % end % % if strcmp(mode,'exp_test') % f = buf(1)*exp(-buf(1)^2-buf(2)^2)+1; % end % % if strcmp(mode,'Rosen_Brock') % f = 100*(buf(2)-buf(1).^2).^2+(1-buf(1)).^2; % end % % end
test_function.m檔案
function f = test_function(buf,mode) if strcmp(mode,'Schaffer') f = 0.5-((sin(sqrt(buf(1).^2+buf(2).^2)).^2)-0.5)./(1+0.001.*(buf(1).^2+buf(2).^2)).^2; end % if mode=='Rastrigin' % f = (buf(1).^2-10*cos(2*pi*buf(1))+10) + (buf(2).^2-10*cos(2.*pi.*buf(2))+10); % end if strcmp(mode,'exp_test') f = buf(1)*exp(-buf(1)^2-buf(2)^2)+1; end if strcmp(mode,'Rosen_Brock') f = 100*(buf(2)-buf(1).^2).^2+(1-buf(1)).^2; end end
程式碼執行結果:
1.Schaffer函式(函式是不是很有特點呢)
尋優的值:
2. Rosen_Brock函式(是不是更有特點呢)
尋優的值:
尋優點在x-y平面上的尋優軌跡:
從尋優的結果發現,演算法陷入了局部最優而沒有跳出來,這是一大問題。
2018-12-2 週日。重度霧霾,閒來無事
------------------------------------------------------------end----------------------------------------------------------------------
期待後續有版本更新