本页采用了机器翻译。点击此处可查看英文原文。

使用自定义数据类型的模拟退火进行多处理器调度

打开实时脚本

此示例说明如何使用模拟退火来最小化使用自定义数据类型的函数。这里，自定义模拟退火来解决多处理器调度问题。

多处理器调度问题

多处理器调度问题包括在一组处理器上找到任务的最佳分配。给出了处理器的数量和任务的数量。处理器完成一项任务所需的时间也作为数据提供。每个处理器独立运行，但每个处理器每次只能运行一个作业。我们将所有作业分配给可用的处理器称为“计划”。该问题的目标是确定给定任务集的最短时间表。

首先，我们确定如何将该问题表述为一个自定义数据类型的优化问题，以便 |simulannealbnd| 函数能够解决该问题。我们提出以下方案：首先，让每个任务用 1 到任务总数之间的整数表示。类似地，每个处理器都由 1 到处理器数量之间的整数表示。现在我们可以将给定作业在给定处理器上所需的时间存储在称为“长度”的矩阵中。处理器 "i" 完成任务 "j" 所花费的时间 "t" 将存储在 lengths(i,j) 中。

我们可以用类似的方式来表示时间表。在给定的计划中，行（1 到处理器数量之间的整数）将代表处理器，列（1 到任务数量之间的整数）将代表任务。例如，调度方案 [1 2 3;4 5 0;6 0 0] 表示任务 1、2 和 3 在处理器 1 上执行，任务 4 和 5 在处理器 2 上执行，任务 6 在处理器 3 上执行。

在这里我们定义任务数量、处理器数量和长度数组。各行不同的系数表示不同处理器以不同的速度工作。我们还定义了一个 "sampleSchedule"，它将作为 simulannealbnd 的起点输入。

rng default % for reproducibility
numberOfProcessors = 11;
numberOfTasks = 40;
lengths = [10*rand(1,numberOfTasks);
    7*rand(1,numberOfTasks);
    2*rand(1,numberOfTasks);
    5*rand(1,numberOfTasks);
    3*rand(1,numberOfTasks);
    4*rand(1,numberOfTasks);
    1*rand(1,numberOfTasks);
    6*rand(1,numberOfTasks);
    4*rand(1,numberOfTasks);
    3*rand(1,numberOfTasks);
    1*rand(1,numberOfTasks)];

% Random distribution of task on processors (starting point)
sampleSchedule = zeros(numberOfProcessors,numberOfTasks);
for task = 1:numberOfTasks
    processorID = 1 + floor(rand*(numberOfProcessors));
    index = find(sampleSchedule(processorID,:)==0);
    sampleSchedule(processorID,index(1)) = task;
end

自定义数据类型的模拟退火

默认情况下，模拟退火算法假设决策变量是双精度数据类型来解决优化问题。因此，生成后续点的退火函数假定当前点是 double 类型的向量。然而，如果将 DataType 选项设置为 "custom"，则模拟退火求解器也可以处理涉及任意数据类型的优化问题。您可以使用任何有效的 MATLAB® 数据结构作为决策变量。例如，如果我们希望 simulannealbnd 使用 "sampleSchedule" 作为决策变量，则可以使用整数矩阵指定自定义数据类型。除了将 DataType 选项设置为 "custom" 之外，我们还需要通过 AnnealingFcn 选项提供一个能够生成新点的自定义退火函数。

自定义退火函数

本节介绍如何创建和使用所需的自定义退火函数。如前所述，多处理器调度问题的试验点是处理器（行）和任务（列）的矩阵。多处理器调度问题的自定义退火函数将以作业调度作为输入。退火函数随后将修改该计划并返回一个与温度成比例改变的新计划（如模拟退火的惯例）。自定义退火函数是 mulprocpermute:

function schedule = mulprocpermute(optimValues,problemData)
% MULPROCPERMUTE Moves one random task to a different processor.
% NEWX = MULPROCPERMUTE(optimValues,problemData) generate a point based
% on the current point and the current temperature

% Copyright 2006 The MathWorks, Inc.

schedule = optimValues.x;
% This loop will generate a neighbor of "distance" equal to
% optimValues.temperature.  It does this by generating a neighbor to the
% current schedule, and then generating a neighbor to that neighbor, and so
% on until it has generated enough neighbors.
for i = 1:floor(max(optimValues.temperature))+1
    [nrows ncols] = size(schedule);
    schedule = neighbor(schedule, nrows, ncols);
end

%=====================================================%
function schedule = neighbor(schedule, nrows, ncols)
% NEIGHBOR generates a single neighbor to the given schedule.  It does so
% by moving one random task to a different processor.  The rest of the code
% is to ensure that the format of the schedule remains the same.

row1 = randinteger(1,1,nrows)+1;
col = randinteger(1,1,ncols)+1;
while schedule(row1, col)==0
    row1 = randinteger(1,1,nrows)+1;
    col = randinteger(1,1,ncols)+1;
end
row2 = randinteger(1,1,nrows)+1;
while row1==row2
    row2 = randinteger(1,1,nrows)+1;
end

for j = 1:ncols
    if schedule(row2,j)==0
        schedule(row2,j) = schedule(row1,col);
        break
    end
end

schedule(row1, col) = 0;
for j = col:ncols-1
    schedule(row1,j) = schedule(row1,j+1);
end
schedule(row1,ncols) = 0;
end
%=====================================================%
function out = randinteger(m,n,range)
%RANDINTEGER generate integer random numbers (m-by-n) in range

len_range = size(range,1) * size(range,2);
% If the IRANGE is specified as a scalar.
if len_range < 2
    if range < 0
        range = [range+1, 0];
    elseif range > 0
        range = [0, range-1];
    else
        range = [0, 0];    % Special case of zero range.
    end
end
% Make sure RANGE is ordered properly.
range = sort(range);

% Calculate the range the distance for the random number generator.
distance = range(2) - range(1);
% Generate the random numbers.
r = floor(rand(m, n) * (distance+1));

% Offset the numbers to the specified value.
out = ones(m,n)*range(1);
out = out + r;
end
end

目标函数

我们需要一个解决多处理器调度问题的目标函数。目标函数返回给定计划所需的总时间（即每个处理器执行其任务所花费时间的最大值）。因此，目标函数也需要长度矩阵才能计算总时间。我们将尽力缩短这个总时间。目标函数是 mulprocfitness：

function timeToComplete = mulprocfitness(schedule, lengths)
%MULPROCFITNESS determines the "fitness" of the given schedule.
%  In other words, it tells us how long the given schedule will take using the
%  knowledge given by "lengths"

%   Copyright 2006 The MathWorks, Inc.

[nrows ncols] = size(schedule);
timeToComplete = zeros(1,nrows);
for i = 1:nrows
    timeToComplete(i) = 0;
    for j = 1:ncols
        if schedule(i,j)~=0
            timeToComplete(i) = timeToComplete(i)+lengths(i,schedule(i,j));
        else
            break
        end
    end
end
timeToComplete = max(timeToComplete);
end

simulannealbnd 仅使用一个参量 x 调用目标函数，但该目标函数有两个参量：x 和 "lengths"。使用一个匿名函数来捕获额外参量（长度矩阵）的值。创建一个函数句柄 ObjectiveFcn，指向一个匿名函数，该函数接受一个输入 x，但会使用 x 和 lengths 调用 mulprocfitness。在创建函数句柄 FitnessFcn 时，变量 lengths 已有值，因此这些值会被匿名函数捕获。

lengths 之前已定义，因此已存在于工作区中。

fitnessfcn = @(x) mulprocfitness(x,lengths);

绘图函数

添加一个自定义绘图函数，用于绘制各处理器上任务所耗费的时间。每个条形代表一个处理器，每个条形的不同颜色块代表不同的任务。

function stop = mulprocplot(~,optimvalues,flag,lengths)
%MULPROCPLOT PlotFcn used for SAMULTIPROCESSORDEMO
%   STOP = MULPROCPLOT(OPTIONS,OPTIMVALUES,FLAG) where OPTIMVALUES is a
%   structure with the following fields:
%              x: current point
%           fval: function value at x
%          bestx: best point found so far
%       bestfval: function value at bestx
%    temperature: current temperature
%      iteration: current iteration
%      funccount: number of function evaluations
%             t0: start time
%              k: annealing parameter "k"
%
%   FLAG: Current state in which PlotFcn is called. Possible values are:
%           init: initialization state
%           iter: iteration state
%           done: final state
%
%   STOP: A boolean to stop the algorithm.
%

%   Copyright 2006-2025 The MathWorks, Inc.

persistent thisTitle %#ok

stop = false;
switch flag
    case "init"
        set(gca,xlimmode="manual",zlimmode="manual", ...
            alimmode="manual")
        titleStr = sprintf("Current Point - Iteration %d", optimvalues.iteration);
        thisTitle = title(titleStr,interp="none");
        toplot = i_generatePlotData(optimvalues, lengths);
        ylabel("Time",interp="none");
        bar(toplot,"stacked",edgecolor="none");
        Xlength = size(toplot,1);        
        set(gca,xlim=[0,1 + Xlength])
    case "iter"
        if ~rem(optimvalues.iteration, 100)
            toplot = i_generatePlotData(optimvalues, lengths);
            bar(toplot,"stacked",edgecolor="none");
            titleStr = sprintf("Current Point - Iteration %d", optimvalues.iteration);
            thisTitle = title(titleStr,interp="none");            
        end
end

    function toplot = i_generatePlotData(optimvalues, lengths)

        schedule = optimvalues.x;
        nrows = size(schedule,1);
        % Remove zero columns (all processes are idle)
        maxlen = 0;
        for i = 1:nrows
            if max(nnz(schedule(i,:))) > maxlen
                maxlen = max(nnz(schedule(i,:)));
            end
        end
        schedule = schedule(:,1:maxlen);

        toplot = zeros(size(schedule));
        [nrows, ncols] = size(schedule);
        for i = 1:nrows
            for j = 1:ncols
                if schedule(i,j)==0 % idle process
                    toplot(i,j) = 0;
                else
                    toplot(i,j) = lengths(i,schedule(i,j));
                end
            end
        end
    end
end

第二个自定义绘图函数

在模拟退火算法中，当前的调度方案未必是迄今为止找到的最佳调度方案。因此，请创建第二个自定义绘图函数，用于显示迄今为止发现的最佳排程方案。

function stop = mulprocplotbest(~,optimvalues,flag,lengths)
%MULPROCPLOTBEST PlotFcn used for SAMULTIPROCESSORDEMO
%   STOP = MULPROCPLOTBEST(OPTIONS,OPTIMVALUES,FLAG) where OPTIMVALUES is a
%   structure with the following fields:
%              x: current point
%           fval: function value at x
%          bestx: best point found so far
%       bestfval: function value at bestx
%    temperature: current temperature
%      iteration: current iteration
%      funccount: number of function evaluations
%             t0: start time
%              k: annealing parameter "k"
%
%   FLAG: Current state in which PlotFcn is called. Possible values are:
%           init: initialization state
%           iter: iteration state
%           done: final state
%
%   STOP: A boolean to stop the algorithm.
%

%   Copyright 2006-2025 The MathWorks, Inc.

persistent thisTitle %#ok

stop = false;
switch flag
    case "init"
        set(gca,xlimmode="manual",zlimmode="manual", ...
            alimmode="manual")
        titleStr = sprintf("Current Point - Iteration %d", optimvalues.iteration);
        thisTitle = title(titleStr,interp="none");
        toplot = i_generatePlotData(optimvalues, lengths);
        Xlength = size(toplot,1);
        ylabel("Time",interp="none");
        bar(toplot,"stacked",edgecolor="none");
        set(gca,xlim=[0,1 + Xlength])
    case "iter"
        if ~rem(optimvalues.iteration, 100)
            toplot = i_generatePlotData(optimvalues, lengths);
            bar(toplot,"stacked",edgecolor="none");
            titleStr = sprintf("Best Point - Iteration %d", optimvalues.iteration);
            thisTitle = title(titleStr,interp="none");            
        end

end

    function toplot = i_generatePlotData(optimvalues, lengths)

        schedule = optimvalues.bestx;
        nrows = size(schedule,1);
        % Remove zero columns (all processes are idle)
        maxlen = 0;
        for i = 1:nrows
            if max(nnz(schedule(i,:))) > maxlen
                maxlen = max(nnz(schedule(i,:)));
            end
        end
        schedule = schedule(:,1:maxlen);

        toplot = zeros(size(schedule));
        [nrows, ncols] = size(schedule);
        for i = 1:nrows
            for j = 1:ncols
                if schedule(i,j)==0
                    toplot(i,j) = 0;
                else
                    toplot(i,j) = lengths(i,schedule(i,j));
                end
            end
        end
    end
end

模拟退火选项设置

选择自定义退火和绘图函数，并修改部分默认选项。将 ReannealInterval 设为 800，因为当求解器开始取得大量局部进展时，ReannealInterval 的较低值似乎会导致温度升高。将 StallIterLimit 值减小到 800，因为默认值会导致求解器运行过慢。最后，将 DataType 设置为 "custom"。

options = optimoptions(@simulannealbnd,DataType="custom", ...
    AnnealingFcn= @mulprocpermute,MaxStallIterations=800,ReannealInterval=800, ...
    PlotFcn={{@mulprocplot, lengths},{@mulprocplotbest, lengths},@saplotf,@saplotbestf});

调用求解器并显示解

问题和选项已全部列出。调用 simulannealbnd 以求解问题。

schedule = simulannealbnd(fitnessfcn,sampleSchedule,[],[],options);

Figure Simulated Annealing contains 4 axes objects. Axes object 1 with title Current Point - Iteration 5400, ylabel Time contains 7 objects of type bar. Axes object 2 with title Best Point - Iteration 5400, ylabel Time contains 8 objects of type bar. Axes object 3 with title Current Function Value: 1.38832, xlabel Iteration, ylabel Function value contains an object of type scatter. Axes object 4 with title Best Function Value: 1.38832, xlabel Iteration, ylabel Function value contains an object of type scatter.

simulannealbnd stopped because the change in best function value is less than options.FunctionTolerance.

从解中移除零列，即所有进程处于空闲状态的列。1

maxlen = 0;
for i = 1:size(schedule,1)
    if max(nnz(schedule(i,:)))>maxlen
        maxlen = max(nnz(schedule(i,:)));
    end
end

显示生成的日程安排。

schedule = schedule(:,1:maxlen)
schedule = 11×8

    22    34    32     0     0     0     0     0
     5     0     0     0     0     0     0     0
    19     6    12    11    39    35     0     0
     7    20     0     0     0     0     0     0
    30    15    10     3     0     0     0     0
    18    28     0     0     0     0     0     0
    31    33    29     4    21     9    25    40
    24    26    14     0     0     0     0     0
    13    16    23     0     0     0     0     0
    38    36     1     0     0     0     0     0
     8    27    37    17     2     0     0     0
      ⋮

另请参阅

主题