智能算法作业

模拟退火代码

clear
%求解电路
syms i1 i2 i3 r1 r2 u1 u2 u3
eqn1 = 2 * (i1 - i3) + r1 * (i1 - i2) + 4 * i1 == u1 - u2;
eqn2 = 3 * i3 + 2 * (i3 - i2) + 2 * (i3 - i1) == 0;
eqn3 = 2 * (i2 - i3) + i2 * r2 + r1 * (i2 - i1) == u2 - u3;
eqns = [eqn1, eqn2, eqn3];
[I1, I2, I3] = solve(eqns, [i1,i2,i3]);
P(u1, u2, u3, r1, r2) = r1 * (I1 - I2) ^ 2 + r2 * I2 ^ 2;

%设置退火参数
alpha = 0.95;
beta = 0.01;
T0 = 180;
Tf = 1;
T = T0;
Markov_length = 100;
coolingTime = ceil((log(Tf) - log(T0))/log(alpha));

%设置初始值
sol_new = [rand()*16, rand()*16, rand()*16, rand()*3+2, rand()*3+2];
sol_current = sol_new;
sol_best = sol_new;
E_new = 0;
E_current = 0;
E_best = 0;
E_bestHistory = zeros(1, coolingTime);
E_currentHistory = zeros(1, coolingTime);

%循环 
loopTimes = 0;
while T >= Tf
    loopTimes = loopTimes + 1;
    for r = 1:Markov_length
       sol_new = changePlan(sol_current);
       E_new = double(P(sol_new(1), sol_new(2), sol_new(3), ...
           sol_new(4),sol_new(5)));
       if E_new > E_current
            E_current = E_new;
            sol_current = sol_new;
            if E_new > E_best
                E_best = E_new;
                sol_best = sol_new;
            end
        else
            if rand < exp((E_new - E_current) ./T ./ beta)
                E_current = E_new;
                sol_current = sol_new;
            end
        end
    end
    T = T.*alpha;
    disp(T)
    E_bestHistory(1, loopTimes) = E_best;
    E_currentHistory(1, loopTimes) = E_current;
end

figure()
plot(E_bestHistory, 'r-')
hold on
plot(E_currentHistory, 'b-')
fh = FigureHandle2D();
fh.plotDecoration('降温次数', '目标函数', '退火效果图', ...
    {'历史最优解', '当前最优解'})
fh.setFontSize(16)

function newA = changePlan(A)
%   扰动函数
    paramLength = length(A);
    newA = A;
    %   扰动 1—3 次
    changeTime = randi(3);
    for i = 1:changeTime
        %   改变的元素位置随机
        changeIndex = randi(paramLength);
        %   扰动范围 -1 ~ 1
        newA(changeIndex)= newA(changeIndex) + rand() * 2 - 1;
    end
    %   使元素满足条件
    %   电压
    if newA(1) + newA(2) + newA(3) > 16
        Usum = newA(1) + newA(2) + newA(3);
        for j = 1:3
            newA(j) = newA(j) / Usum * 16;
        end
    end
    [~, UmaxIndex] = max(newA(1:3));
    for i = 1:3
        if newA(i) < 1
            newA(UmaxIndex) =  newA(UmaxIndex) - (1 - newA(i));
            newA(i) = 1;
        end
    end
    %   电阻
    for i = 4:5
       if newA(i) < 2 || newA(i) > 5
            newA(i) = A(i);
       end
    end
end

运算结果

U1	U2	U3	R1	R2	P
1	14	1	2.0067	2.0005	24.6938

遗传算法

clear
%   先解电路
syms i1 i2 i3 r1 r2 u1 u2 u3
eqn1 = 2 * (i1 - i3) + r1 * (i1 - i2) + 4 * i1 == u1 - u2;
eqn2 = 3 * i3 + 2 * (i3 - i2) + 2 * (i3 - i1) == 0;
eqn3 = 2 * (i2 - i3) + i2 * r2 + r1 * (i2 - i1) == u2 - u3;
eqns = [eqn1, eqn2, eqn3];
[I1, I2, I3] = solve(eqns, [i1, i2, i3]);
P(u1, u2, u3, r1, r2) = r1 * (I1 - I2) ^ 2 + r2 * I2 ^ 2;

%   设置遗传变异参数
popSize = 50;
binaryLength = 20;
mutationPosibility = 0.007;
crossPosibility = 0.5;
mutationTimes = 100;

%   初始化种群、初始化记录数组
pop = round(rand(popSize, binaryLength * 4));
resultHistory = zeros(popSize, mutationTimes + 1);
bestResult = 0;
bestParameter = zeros(1, 5);

%   开始进化
for index = 1:mutationTimes
    disp(index);
    %   交叉、变异
    pop = crossover(pop, crossPosibility);
    pop = mutate(pop, mutationPosibility);
    %   求 pop 对应参数和结果
    decodePop = decodechrom(pop,binaryLength);
    paramArray = getParam(decodePop, binaryLength);
    [result, fitValue] = getFitValue(paramArray, P);
    %   存储记录
    resultHistory(:, index) = result;
    if max(result) > bestResult
        [bestResult, bestIndex] = max(result);
        bestParameter = paramArray(bestIndex,:);
    end
    %   选择
    pop = select(pop, fitValue);
end
%   计算最后一次选择后的结果
decodePop = decodechrom(pop,binaryLength);
paramArray = getParam(decodePop, binaryLength);
[result, fitValue] = getFitValue(paramArray, P);
resultHistory(:, index + 1) = result;

%   作图
figure()
bestHistory = max(resultHistory);
avgHistory = sum(resultHistory) / size(resultHistory, 1);
plot(avgHistory, 'b-');
hold on
plot(bestHistory, 'r-');

function paramArray = getParam(decodePop, binaryLength)
%   获取参数，[u1, u2, u3, r1, r2]
%   decodePop 为 [u1, u2, r1, r2]
    paramArray = zeros(size(decodePop, 1), 5);
    maxNum = 2 ^ binaryLength - 1;
    for index = 1 : size(decodePop, 1)
       u1 = decodePop(index, 1) / maxNum * 14 + 1;
       u2 = decodePop(index, 2) / maxNum * 14 + 1;
       r1 = decodePop(index, 3) / maxNum * 3 + 2;
       r2 = decodePop(index, 4) / maxNum * 3 + 2;
       if u1 + u2 > 15
           u1 = 15 * u1 / (u1 + u2);
           u2 = 15 - u1;
           u3 = 1;
       else
           u3 = 16 - u1 - u2;
       end
       paramArray(index, :) = [u1, u2, u3, r1, r2];
    end
end

function [result, fitValue] = getFitValue(paramArray, P)
%   获取功率和适应度
    result = zeros(size(paramArray, 1), 1);
    for index = 1 : size(paramArray, 1)
        result(index, 1) = P(paramArray(index, 1), paramArray(index, 2) ...
            ,paramArray(index, 3), paramArray(index, 4) ...
            ,paramArray(index, 5));
    end
    fitValue = result;
end

%   以下函数为遗传算法函数，通用
function pop2 = decodebinary(pop)
%   将二进制数转换为十进制数
    popColumn = size(pop,2);
    pop2 = zeros(size(pop,1),1);
    for i = 1:popColumn
        pop2 = pop2 + pop(:,i) * 2^(popColumn - i);
    end
end

function pop2 = decodechrom(pop, binaryLength)
%   将二进制数组转换为十进制数组
    popColumn = size(pop,2);
    decimalLength = popColumn/binaryLength;
    if decimalLength ~= round(decimalLength)
        error('数组列数错误');
    end
    pop2 = zeros(size(pop,1),decimalLength);
    for i = 1:decimalLength
        beginNumber = binaryLength * (i-1) + 1;
        endNumber = binaryLength * i;
        pop2(:,i) = decodebinary(pop(:,beginNumber:endNumber));
    end
end

function pop2 = mutate(pop, mutationPosibility)
%   变异
    number = size(pop,1) * size(pop,2);
    mutationNumber = round(mutationPosibility * number);
    pop2 = pop;
    for i = 1:mutationNumber
        indexRow = randi(size(pop,1));
        indexColumn = randi(size(pop,2));
        pop2(indexRow,indexColumn) = double(~pop2(indexRow,indexColumn));
    end
end

function pop2 = crossover(pop, crossPosibility)
%   交叉
    pop2 = pop;
    for i = 1:size(pop,1) - 1
        if rand() < crossPosibility
            index = randi(size(pop,2));
            pop2(i,index:end) = pop(i + 1,index:end);
            pop2(i + 1,index:end) = pop(i,index:end);
        end
    end
end

function pop2 = select(pop, fitValue)
%   选择
    pop2 = zeros(size(pop));
    fitValue = cumsum(fitValue) / sum(fitValue);
    for i = 1:size(pop2,1)
        index = find(fitValue>=rand(),1);
        pop2(i,:) = pop(index,:);
    end
end

运算结果

U1	U2	U3	R1	R2	P
1.0159	13.9821	1.0020	2.1037	2.0398	24.2664

根据两个算法，可以看出找到的解都是接近的，由于这个问题自变量很少，而且存在智能算法的精度误差，对各参数取整可以得到更好的解：

U1	U2	U3	R1	R2	P
1	14	1	2	2	24.7116

---

本文章使用limfx的vsocde插件快速发布