function ret = u(nk, nm, np, lam)
  ret = 1 / (1/5 - lam/nk) + 1 / (1/23 - lam/(3 * np)) + 1 / (1/13 - lam/(3 * nm));
endfunction

function ret = ct(nk, nm, np, lam)
  ret = u(nk, nm, np, lam) * (3 * 15 * np + 3 * 6 * nm + 2 * nk);
endfunction

lambda = 3/23;
b = [0.2, 0.4, 0.3, 0.1];
a = [0.1, 0.4, 0.7, 0.9];

ma = []; mb = []; mc = [];
m1 = m2 = m3 = m4 = [];
for nk = 1:93
    for nm = 1:11
        for np = 1:5
            tbl = table(nk, nm, np);
            if (tbl(3) <= 250)
                acc = 0;
                for i = 1:4
                    acc = acc + b(i) * ct(nk, nm, np, lambda * a(i));
                endfor
                ma = [ms; [nk, nm, np, acc]];
                mb = [mb; [nk, nm, np, ct(nk, nm, np, 0.48 * lambda)]];
                mc = [mc; [nk, nm, np, ct(nk, nm, np, lambda - 0.01)]];
                m1 = [m1; [nk, nm, np, ct(nk, nm, np, lambda * a(1))]];
                m2 = [m2; [nk, nm, np, ct(nk, nm, np, lambda * a(2))]];
                m3 = [m3; [nk, nm, np, ct(nk, nm, np, lambda * a(3))]];
                m4 = [m4; [nk, nm, np, ct(nk, nm, np, lambda * a(4))]];
            endif
        endfor
    endfor
endfor
[min_vals idc] = min(ma);
ma(idc(4), :)
[min_vals idc] = min(mb);
mb(idc(4), :)
[min_vals idc] = min(mc);
mc(idc(4), :)
