的最小值，其中每個變數的取值區間都是[-1, +1]。

首先，確定適應度的函式：

function y = my_fitness(population)
% population是隨機數[0,1]矩陣，下面的操作改變範圍為[-1,1]
population = 2 * (population - 0.5); 
y = sum(population.^2, 2); % 行的平方和

function [best_fitness, elite, generation, last_generation] = my_ga( ...
    number_of_variables, ...    % 求解問題的引數個數
    fitness_function, ...       % 自定義適應度函式名
    population_size, ...        % 種群規模（每一代個體數目）
    parent_number, ...          % 每一代中保持不變的數目（除了變異）
    mutation_rate, ...          % 變異概率
    maximal_generation, ...     % 最大演化代數
    minimal_cost ...            % 最小目標值（函式值越小，則適應度越高）
)

% 累加概率
% 假設 parent_number = 10
% 分子 parent_number:-1:1 用於生成一個數列
% 分母 sum(parent_number:-1:1) 是一個求和結果（一個數）
% 分子 10     9     8     7     6     5     4     3     2     1
% 分母 10+9+...+1=55
% 相除 0.1818    0.1636    0.1455    0.1273    0.1091    0.0909    0.0727    0.0545    0.0364    0.0182
% 累加(cumsum) 0.1818    0.3455    0.4909    0.6182    0.7273    0.8182    0.8909    0.9455    0.9818    1.0000
% 運算結果可以看出: 累加概率函式是一個從0到1增長得越來越慢的函式
cumulative_probabilities = cumsum((parent_number:-1:1) / sum(parent_number:-1:1)); % 1個長度為parent_number的數列
best_fitness = ones(maximal_generation, 1);% 記錄每一代最佳適應度,先初始化為1
elite = zeros(maximal_generation, number_of_variables);% 記錄每一代的最優解，初始化為0

% 子女數量=種群數量 - 父母數量（父母即每一代中不發生改變的個體）
child_number = population_size - parent_number; % 每一代子女的數目

% population_size 對應矩陣的行，每一行表示1個個體，行數=個體數（種群數量）
% number_of_variables 對應矩陣的列，列數=引數個數（個體特徵由這些引數表示）
population = rand(population_size, number_of_variables);% 初始化種群
last_generation = 0; % 記錄跳出迴圈時的代數

for generation = 1 : maximal_generation % 演化迴圈開始
    
    % feval把資料帶入到一個定義好的函式控制代碼中計算
    % 把population矩陣帶入fitness_function函式計算
    cost = feval(fitness_function, population); % 計算所有個體的適應度（population_size*1的矩陣）

    % index記錄排序後每個值原來的行數
    [cost, index] = sort(cost); % 將適應度函式值從小到大排序

    % index(1:parent_number) 
    % 前parent_number個cost較小的個體在種群population中的行數
    % 選出這部分(parent_number個)個體作為父母，其實parent_number對應交叉概率
    population = population(index(1:parent_number), :); % 先保留一部分較優的個體
    % population矩陣是不斷變化的
    % cost在經過前面的sort排序後，矩陣已經改變為升序的
    % cost(1)即為本代的最佳適應度
    best_fitness(generation) = cost(1); % 記錄本代的最佳適應度

    % population矩陣第一行為本代的最優解（精英）
    elite(generation, :) = population(1, :); % 記錄本代的最優解（精英）

    % 若本代的最優解已足夠好，則停止演化
    if best_fitness(generation) < minimal_cost; 
        last_generation = generation;
        break; 
    end
    
    % 交叉變異產生新的種群，染色體交叉開始
    for child = 1:2:child_number % 步長為2是因為每一次交叉會產生2個孩子
        
        % cumulative_probabilities 長度為 parent_number
        % 從中隨機選擇2個父母出來  (child+parent_number)
        mother = find(cumulative_probabilities > rand, 1); % 選擇一個較優秀的母親
        father = find(cumulative_probabilities > rand, 1); % 選擇一個較優秀的父親
        
        % ceil：向上取整
        % rand 生成一個隨機數
        % 即隨機選擇了一列，這一列的值交換
        crossover_point = ceil(rand*number_of_variables); % 隨機地確定一個染色體交叉點
        
        % 假如crossover_point=3, number_of_variables=5
        % mask1 = 1     1     1     0     0
        % mask2 = 0     0     0     1     1
        mask1 = [ones(1, crossover_point), zeros(1, number_of_variables - crossover_point)];
        mask2 = not(mask1);
        
        % 獲取分開的4段染色體
        mother_1 = mask1 .* population(mother, :); % 母親染色體的前部分
        mother_2 = mask2 .* population(mother, :); % 母親染色體的後部分
        
        father_1 = mask1 .* population(father, :); % 父親染色體的前部分
        father_2 = mask2 .* population(father, :); % 父親染色體的後部分
        
        % 得到下一代
        population(parent_number + child, :) = mother_1 + father_2; % 一個孩子
        population(parent_number+child+1, :) = mother_2 + father_1; % 另一個孩子
        
    end % 染色體交叉結束
    
    
    % 染色體變異開始，變異種群
    mutation_population = population(2:population_size, :); % 精英不參與變異，從2開始    
    number_of_elements = (population_size - 1) * number_of_variables; % 全部基因數目
    number_of_mutations = ceil(number_of_elements * mutation_rate); % 變異的基因數目（基因總數*變異率）
    
    % rand(1, number_of_mutations) 生成number_of_mutations個隨機數(範圍0-1)組成的矩陣(1*number_of_mutations)
    % 數乘後，矩陣每個元素表示發生改變的基因的位置（元素在矩陣中的一維座標）
    mutation_points = ceil(number_of_elements * rand(1, number_of_mutations)); % 確定要變異的基因
    
    % 被選中的基因都被一個隨機數替代，完成變異
    mutation_population(mutation_points) = rand(1, number_of_mutations); % 對選中的基因進行變異操作
    population(2:population_size, :) = mutation_population; % 發生變異之後的種群
    % 染色體變異結束
end % 演化迴圈結束
 

clc,clear all,close all;

% 呼叫 my_ga 進行計算
% 求解問題的引數個數         10
% 自定義適應度函式名         my_fitness
% 種群規模                  100
% 每一代中保持不變的數目     50 (即交叉率0.5)
% 變異概率                  0.1 (1/10的個體發生變異)
% 最大演化代數              10000 10000代
% 最小目標值                1.0e-6 個體適應度函式值 < 0.000001結束
[best_fitness, elite, generation, last_generation] = my_ga(10, 'my_fitness', 100, 50, 0.1, 10000, 1.0e-6);

% 輸出後10行
disp(last_generation); 
i_begin = last_generation - 9;
disp(best_fitness(i_begin:last_generation,:));
% disp(best_fitness(9990:10000,:));
% disp(elite(9990:10000,:))
% 不要寫成這樣，因為GA常常在中間就跳出迴圈了

% 將elite值轉化為問題範圍內
my_elite = elite(i_begin:last_generation,:);
my_elite = 2 * (my_elite - 0.5);
disp(my_elite);

% 最佳適應度的演化情況
figure
loglog(1:generation, best_fitness(1:generation), 'linewidth',2)
xlabel('Generation','fontsize',15);
ylabel('Best Fitness','fontsize',15);
set(gca,'fontsize',15,'ticklength',get(gca,'ticklength')*2);

% 最優解的演化情況
figure
semilogx(1 : generation, 2 * elite(1 : generation, :) - 1)
xlabel('Generation','fontsize',15);
ylabel('Best Solution','fontsize',15);
set(gca,'fontsize',15,'ticklength',get(gca,'ticklength')*2);