6

First of all, i'd like to come clean and say that the following question is for school so don't be too harsh on me :)

I'm having a bit of a problem modelling an optimization problem in matlab using a recursive algorithm (which is a requirement).

The problem's definition is:

Decide the quantity of fish to catch each year considering a time window of 10 years knowing there are presently 10000 fishes in the lake, year 1, the growing rate of fish is the number of fishes present in the lake at the beginning of each year + 20%.

Let x be the quantity of fish to catch, $5 the price of each fish and the cost of catching fish: 

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

decide the number of fish to catch each year, for 10 years, in order to maximize the profit.

Future gains are depreciated by a factor of 0.2/year, which means that earning $1 in year 1 is the same as $0.8 in year 2 and so on.

I currently have defined the following objective function:

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)}

How would i go about implementing this in matlab?

This is what i have so far:

Main file

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

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;

This only computes for the case where the year is 10 and b = 10 (value in thousands). This is the same poblem describes as "Discounted Profits Problem" in book

Any help you can give me will be greatly appreciated.

EDIT 1: I'm really stuck here guys. If you could help me implement this in say Java i would try and port it to Matlab.

EDIT 2: I edited the code to the most recent version. Now i'm getting

"Maximum recursion limit of 500 reached."

Can you help me?

EDIT 3: I managed to get it working but it only returns 0.

EDIT 4: Code updated. Now i'm getting

Attempted to access prop(2,0); index must be a positive integer or logical.

Error in Main (line 66) kq=prop(q,y)

Improve this question
8
  • 1
    I can understand if you are unable to implement a dynamic programming algorithm, this is a bit tricky. But there are some parts of the solution you can definitely provide. Provide code for c(x,b). Provide code for the gain, this should be something like g(x,y) with y=year. Commented Nov 3, 2013 at 14:33
  • 1
    A small comment: in pesca_rec you keep calling pesca_rec(Year+1,...) that doesn't seem to terminate - probably will give you errors. Commented Nov 4, 2013 at 14:28
  • 1
    Also consider making MaxProfit not global cause pesca_rec you keep calling pesca_rec(Year+1,...) and both will write to MaxProfit - so I wonder if it contains the value you intend it to have. Commented Nov 4, 2013 at 14:32
  • Please make sure to always run code with dbstop if error the error that you get for edit4 is trivial if you do this and inspect the relevant variables. Commented Nov 5, 2013 at 17:33
  • 1
    Are you sure that you have your cost function right? Usually a semi-real-world cost function like this would be continuous, or nearly so, and show economies of scale. Commented Nov 8, 2013 at 16:35

2 Answers 2

Reset to default
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

After reading my Solution again I have some Ideas to improve the performance:

  • Instead of indexing vcosts with a vector (possibleToFish), directly use fishes to index.
  • Preallocate options / crate in one step

For 10000 it runs in acceptable time (about 5 min), for bigger data I would recommend to update.

Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

0

Relatively clean and classic implementation of a solution using dynamic programming, that should give the guaranteed optimal solution (assuming no bugs):

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 with this is that it is pretty slow, I guess that the running time is something like O(nyear * (nstart*grow_rate^(nyear-1))^2). This is polynomial time (?), except for the exponential grow_rate. (Note that, due to the caching, this is still better than a brute-force solution, which would be O((nstart*grow_rate^(nyear-1))^nyear), which is exponential in nyear.) The quadratic dependence on nstart makes it a bit too slow to run for nstart = 10000 (probably around a day), but for nstart = 200 it still runs in acceptable time. Some 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 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

So the depreciated rate is so high that the optimal solution is to empty the complete lake in the first year. Reducing this depreciated rate a bit (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 growing rate is higher than the devaluation rate, so the optimal solution is to let the fish procreate as long as possible, and empty out the lake in the last year. Interestingly, the solution proposes to catch a handful of fish in the intermediate years. I assume that these extra fish would not make a change in the rounding in nleft = floor(grow_rate * nleft), so you gain a tiny bit by fishing them out earlier. Playing a bit with dep_rate, it seems a pretty critical all-or-nothing situation: if the grow_rate is high enough, you fish everything in the last year. If it is too low, you fish everything in the first year.

Improve this answer

4 Comments

Your initial value of fishes is 200 instead of 10000, this makes the cost function linar and changes the results.
I know, the results are only optimal under those conditions. But I don't expect that this will change the all-or-nothing effect a lot, unless the cost function is strongly non-linear.
It changes, the optimal solution seems to fish 5000,5000,4080. At least when I run my solution.
Ah yes, I guess this is because there is a large step in the cost at 5000, so it is cheaper to fish 5000 at a time. This is is not really realistic, as Pursuit noted. I could try to introduce a step around 100 or so in my cost function, but I don't have time now.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.