动态规划 - 递归实现

Question

首先，我想坦白地说，以下问题是为学校而提出的，请不要对我太苛刻 :)

我在使用 MATLAB 中的递归算法（这是一个要求）建模优化问题时遇到了一些问题。

问题的定义是：

考虑一个为期10年的时间窗口，决定每年捕捞多少鱼，已知湖中目前有10000条鱼，第1年中鱼的增长率为每年开始时湖中鱼的数量的20%。

设 x 为捕捉的鱼的数量，每条鱼的价格为 $5，捕鱼的成本为：

0.4x + 100 if x is <= 5000; 
0.3x + 5000 if  5000 <= x <= 10000; 
0.2x + 10000 if x > 10000; 

为了最大化利润，在未来10年中决定每年捕捞的鱼数量。

未来收益将按照每年0.2的因子递减，这意味着在第一年赚取1美元与第二年赚取0.8美元相同，以此类推。

目前我已经定义了以下目标函数：

x -> quantity of fish to catch
b-> quantity of fish availavable in the beginning of year i
c(x,b) -> cost of catching x fish with b fishes available

f_i(b) = max {(5x - c(x,b)) + 0.8 * f_i+1((b - x) * 1.2)}

我如何在Matlab中实现这个功能？

目前，我已经实现了以下内容：

主文件

clear;

global M Kdep Cost RecursiveProfit ValorF prop

Kdep=[10; 20; 30; 40; 50; 60; 70; 80; 90; 100]; %max nr of fish in the lake at the beginning of each year, 10 years, in thousands. Growth factor = 20%

M=1000;

%Cost and Profit of selling i fishes given that there are j at the beginning of the year
for i = 1:50
    for j = 1:11
        Cost(i,j) = 0.2 * i + 10;
        RecursiveProfit(i,j) = 5 * i - Cost(i, j);
    end
end


for i = 1:10
    for j = 1:10
        Cost(i,j) = 0.3 * i + 5;
        RecursiveProfit(i,j) = 5 * i - Cost(i, j);
    end
end

for i = 1:5
    for j = 1:5
        Cost(i,j) = 0.4 * i + 0.1;
        RecursiveProfit(i,j) = 5 * i - Cost(i, j);
    end
end

%prop = 1 : 10;

ValorF = -M * ones(10, 50);


for a = 1:5
    ValorF(10, a) = 5 * a - (0.4 * a + 1); %On Year 10, if there are <= a thousand fishes in the lake ...
    prop(10, a) = a;
end

for b = 6:10
    ValorF(10, b) = 5 * b - (0.3 * b + 5); %On Year 10, if there are 6 <= a <= 10  thousand fishes in the lake ...
    prop(10, b) = b;
end

for c = 10:41
    ValorF(10, c) = 5 * c - (0.2 * c + 10); 
    prop(10, c) = c;
end

MaxProfit = RecursiveProfit(1, 10)

k1 = prop(1,10)

kant=k1;

y = 6 - Cost(kant,10);

for q=2:10
    if(kant == 0)
    kant = kant + 1;
end
    kq=prop(q,y)
    kant=kq;
    y = y - Cost(kant,q);
end %for i

函数

function y=RecursiveProfit(j,x)
global M Kdep Cost Prof ValorF prop

y=ValorF(j,x);

if y~= -M
    return
end %if

auxMax=-M;
decision=0;

for k=1:Kdep(j)
    if Prof(k,j) <= x-k
        aux=Prof(k,j)+RecursiveProfit(j+1, (x - k));
            if auxMax < aux 
                auxMax=aux;
                decision=k;
            end %if aux
        else break
    end %if Cost   

end %for k

ValorF(j,x)=auxMax;
prop(j,x)=decision;
y=auxMax;

这仅适用于年份为10且b = 10（以千为单位）的情况。这与该书中描述的“折扣利润问题”相同。

任何你能给我的帮助都将不胜感激。

编辑1：我真的被卡住了。如果您可以帮助我在Java中实现它，我会尝试将其移植到Matlab中。

编辑2：我已经更新了代码，现在我得到了“最大递归深度500”的错误。你能帮我吗？

编辑3：我设法让它工作了，但它只返回0。

编辑4：代码已更新。现在我得到了“尝试访问prop（2,0）;索引必须是正整数或逻辑。”的错误。

- DaveQuinn

---

我能理解如果你无法实现动态规划算法，这有点棘手。但是你肯定可以提供解决方案的某些部分。提供c(x,b)的代码。提供增益的代码，这应该类似于g(x,y)，其中y=年份。 - Daniel

一个小注释：在 pesca_rec 中，你不断地调用 pesca_rec(Year+1,...)，这似乎没有终止 - 可能会导致错误。 - bdecaf

还要考虑将MaxProfit不设为全局变量，因为你一直在调用pesca_rec(Year+1,...)，两者都会写入MaxProfit - 所以我想知道它是否包含你想要的值。 - bdecaf

捕获5000条鱼的成本是多少？2100还是6500？你可能需要将一个<=改为<。 - Bas Swinckels

你确定你的成本函数是正确的吗？通常这样的半真实世界的成本函数应该是连续的，或者几乎连续的，并显示规模经济。 - Pursuit

Answer 1

function gofishing(numoffishes,years)

growFactor=1.2;
%index shift, index 1 : 0 fishes
earn{1}=-inf(numoffishes+1,1);
%index shift, index 1 : 0 fishes
earn{1}(numoffishes+1)=0;
%previous: backpointer to find path of solution.
previous{1}=nan;

%index shift, index 1 : 0 fishes
vcosts = zeros(1,ceil(numoffishes*growFactor^years));

for idx=1:numel(vcosts)
    vcosts(idx)=costs(idx-1);
end

for step = 1:years*2
    fprintf('step %d\n',step);
    if mod(step,2)==1;
        %do fish grow step
        earn{step+1}=-inf(floor(numel(earn{step})*1.2)-1,1);
        previous{step+1}=nan(floor(numel(earn{step})*1.2)-1,1);
        for fishes=0:numel(earn{step})-1
            grownfishes=floor(fishes*1.2);
            earn{step+1}(grownfishes+1)=earn{step}(fishes+1);
            previous{step+1}(grownfishes+1)=fishes;
        end
    else
        %do fishing step
        earn{step+1}=-inf(size(earn{step}));
        previous{step+1}=nan(size(earn{step}));
        for fishes=0:numel(earn{step})-1
            if isinf(earn{step}(fishes+1))
                %earn is -inf, nothing to do
                continue;
            end
            possibleToFish=fishes:-1:0;
            %calculate earn for possible amounts to fish
            options=((vrevenue(possibleToFish)-vcosts(possibleToFish+1))*0.8^(step/2-1)+earn{step}(fishes+1))';
            %append -inf for not existing options
            options=[options;-Inf(numel(earn{step+1})-numel(options),1)];
            %found better option:
            better=earn{step+1}<options;
            earn{step+1}(better)=options(better);
            previous{step+1}(better)=fishes;
        end
    end
end
[~,fc]=max(earn{end});
fc=fc-1;
fprintf('ending with %d fishes and a earn of %d\n',fc,earn{end}(fc+1));
for step=(years*2):-1:2
    fc=previous{step}(fc+1);
    fprintf('fish count %d\n',fc');
end
end

function c=costs(x)
if (x<=5000)
    c=0.4*x + 100;
    return
end
if (x <= 10000)
    c=0.3*x + 5000;
    return
end
c=0.2*x + 10000;
return
end
function c=vrevenue(x)
c=5.*x;
end

---

再次阅读我的解决方案后，我有一些可提高性能的想法：

• 不要使用向量（possibleToFish）来索引 vcosts，而是直接使用 fishes 进行索引。
• 一步预分配选项 / 箱

对于10000而言，它可以在可接受的时间内运行（约5分钟），对于更大的数据，建议进行更新。

- Daniel

Answer 2

相对干净和经典的实现方案，使用动态规划，应该提供保证的最优解决方案（假设没有错误）：

function [max_profit, Ncatch] = fish(nstart, nyear, grow_rate, dep_rate)
% return the maximum possible profit max_prof and a vector Ncatch which says
% how many fish to catch for each year

global cache

if isempty(cache) % init_cache
    maxfish = ceil(nstart * grow_rate^nyear);
    % allocate abundant cache space for dynamic programming
    cache.profit = nan(maxfish+1, nyear);
    cache.ncatch = cell(maxfish+1, nyear);
    % function calls are expensive, cache calls to revenue and cost
    cache.netprofit = arrayfun(@(i) revenue(i) - cost(i), 0:maxfish);
    cache.grow_rate = grow_rate;
    cache.dep_rate = dep_rate;
end

if ~isnan(cache.profit(nstart+1, nyear)) % found in cache
    max_profit = cache.profit(nstart+1, nyear);
    Ncatch = cache.ncatch{nstart+1, nyear};

    %assert(cache.grow_rate == grow_rate, 'clear cache first!')
    %assert(cache.dep_rate == dep_rate, 'clear cache first!')
    return 
end

max_profit = -inf;

if nyear == 1 % base case to stop recursion
    % simply get largest profit, future be damned
    [max_profit, imx] = max(cache.netprofit(1:nstart+1));
    Ncatch = [imx - 1];

else % recursive part
    for ncatch = 0:nstart % try all possible actions
        nleft = nstart - ncatch; % catch
        nleft = floor(grow_rate * nleft); % reproduce

        % recursive step, uses optimal profit for 1 year less
        [rec_profit, rec_Ncatch] = fish(nleft, nyear - 1, grow_rate, dep_rate); 

        profit = cache.netprofit(ncatch + 1) + dep_rate * rec_profit;
        if profit > max_profit
            max_profit = profit;
            Ncatch = [ncatch, rec_Ncatch];
        end
    end
end

% store result in cache
cache.profit(nstart+1, nyear) = max_profit;
cache.ncatch{nstart+1, nyear} = Ncatch;
end

function c = cost(x)
    if (x <= 5000)
        c = 0.4 * x + 100;
        return
    end
    if (x <= 10000)
        c = 0.3 * x + 5000;
        return
    end
    c = 0.2 * x + 10000;
end

function r = revenue(x)
    r = 5 .* x;
end

唯一的问题是它运行得相当慢，我猜运行时间大约是 O(nyear * (nstart*grow_rate^(nyear-1))^2)。这是多项式时间（？），除了指数增长率。 （请注意，由于缓存，这仍然比暴力解决方案更好，后者将是 O((nstart*grow_rate^(nyear-1))^nyear)，在 nyear 上呈指数增长。）对 nstart 的二次依赖使其运行速度有点太慢，无法为 nstart = 10000 运行（大约需要一天），但对于 nstart = 200，它仍然可以在可接受的时间内运行。一些快速测试：

clear global % clear the cache

global cache

nyear = 10;
nstart = 200;
grow_rate = 1.2;
dep_rate = 0.80;

tic;
[profit, ncatch] = fish(nstart, nyear, grow_rate, dep_rate);
toc;

nfish = nstart;
for i = 1:nyear
    nnew = floor(grow_rate * (nfish - ncatch(i)));
    prof = cache.netprofit(ncatch(i) + 1);
    dep_prof = prof * (dep_rate)^(i-1);
    fprintf('year %d: start %d, catch %d, left %d, profit %.1f, dep. prof %.1f\n', ...
        i, nfish, ncatch(i), nnew, prof, dep_prof);
    nfish = nnew;
end
fprintf('Total profit: %.1f\n', profit);

带有结果

>> test_fish
Elapsed time is 58.591110 seconds.
year 1: start 200, catch 200, left 0, profit 820.0, dep. prof 820.0
year 2: start 0, catch 0, left 0, profit -100.0, dep. prof -80.0
year 3: start 0, catch 0, left 0, profit -100.0, dep. prof -64.0
year 4: start 0, catch 0, left 0, profit -100.0, dep. prof -51.2
year 5: start 0, catch 0, left 0, profit -100.0, dep. prof -41.0
year 6: start 0, catch 0, left 0, profit -100.0, dep. prof -32.8
year 7: start 0, catch 0, left 0, profit -100.0, dep. prof -26.2
year 8: start 0, catch 0, left 0, profit -100.0, dep. prof -21.0
year 9: start 0, catch 0, left 0, profit -100.0, dep. prof -16.8
year 10: start 0, catch 0, left 0, profit -100.0, dep. prof -13.4
Total profit: 473.6

因此，折旧率非常高，最优解是在第一年完全排空湖泊。将这个折旧率稍微降低一些（dep_rate = 0.84;），情况就会发生逆转：

>> test_fish
Elapsed time is 55.872516 seconds.
year 1: start 200, catch 0, left 240, profit -100.0, dep. prof -100.0
year 2: start 240, catch 0, left 288, profit -100.0, dep. prof -84.0
year 3: start 288, catch 3, left 342, profit -86.2, dep. prof -60.8
year 4: start 342, catch 2, left 408, profit -90.8, dep. prof -53.8
year 5: start 408, catch 3, left 486, profit -86.2, dep. prof -42.9
year 6: start 486, catch 1, left 582, profit -95.4, dep. prof -39.9
year 7: start 582, catch 2, left 696, profit -90.8, dep. prof -31.9
year 8: start 696, catch 1, left 834, profit -95.4, dep. prof -28.2
year 9: start 834, catch 4, left 996, profit -81.6, dep. prof -20.2
year 10: start 996, catch 996, left 0, profit 4481.6, dep. prof 933.1
Total profit: 471.4

在这种情况下，增长率高于贬值率，因此最优解是让鱼尽可能地繁殖，最后一年清空湖泊。有趣的是，该解决方案建议在中间年份捕捉少量的鱼。我认为这些额外的鱼不会改变 nleft = floor(grow_rate * nleft) 中的四舍五入，因此提前捕捞它们可以获得微小的收益。通过调整 dep_rate，似乎是一个非常关键的全有或全无的情况：如果增长率足够高，则在最后一年捕捞所有鱼。如果太低，则在第一年捕捞所有鱼。

- Bas Swinckels

---

你的鱼的初始值是200而不是10000，这使得成本函数变为线性并改变了结果。 - Daniel

我知道，只有在那些条件下才能得到最佳结果。但是除非成本函数非常非线性，否则我不指望这会对全有或全无的效果产生太大影响。 - Bas Swinckels

它会变化，最优解似乎是钓鱼5000,5000,4080。至少在我运行我的解决方案时是这样的。 - Daniel

啊是的，我猜这是因为在5000处成本有一个大跃变，所以一次性捕鱼5000条更便宜。正如Pursuit所指出的那样，这并不真实。我可以尝试在我的成本函数中引入100左右的步长，但我现在没有时间。 - Bas Swinckels