A relatively clean and classic implementation of the solution using dynamic programming, which should provide a guaranteed optimal solution (in the absence of errors):
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
The only problem is that it is rather slow, I think the runtime is something like O(nyear * (nstart*grow_rate^(nyear-1))^2) . This is polynomial time (?), With the exception of exponential grow_rate. (Note that caching is even better than brute force solution, which will be O((nstart*grow_rate^(nyear-1))^nyear) , which exponentially in nyear .) The quadratic dependence on nstart makes it too run nstart = 10000 slowly (maybe about a day), but for nstart = 200 it still works at an acceptable time. Quick test:
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);
with the result
>> 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
Thus, the depreciation rate is so high that the best solution is to clean the entire lake in the first year. The decrease in this depreciated speed is a little ( dep_rate = 0.84; ), the situation turns around:
>> 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
In this case, the growth rate is higher than the devaluation rate, so the best solution is to allow the fish to breed for as long as possible and release the lake last year. Interestingly, the solution suggests catching a few fish in intermediate years. I assume that this extra fish will not change the rounding in nleft = floor(grow_rate * nleft) , so you will gain a little by catching them earlier. Playing a little with dep_rate , this seems like a pretty critical all-or-nothing situation: if grow_rate is tall enough, you catch everything last year. If it is too low, you catch everything in the first year.