1
$\begingroup$

I have the following linear program that I am able to solve in MATLAB. However, I want to move to Mathematica.

For some fixed constants $n$, $\delta$ and $\varepsilon$ and fixed $(n+1)$-dimensional vector $p$, we have the following optimization variables: $(n+1)$-dimensional vectors $y$ and $r$ and a constant $s$. The components of the vectors are labelled by index $k$ which runs from $0$ to $n$. The linear program is

$$ \begin{aligned} \inf \quad & s \\ \text { s.t. } \quad & y_{k}-r_{k}+p_{k} \geq 0, \forall k \\ \quad & y_{k} \geq 0, \forall k \\ \quad& 0 \leq r_{k} \leq s, \forall k \\ & \Sigma_{k=0}^{n}\left(\begin{array}{l} n \\ k \end{array}\right) r_{k}=1, \\ & \Sigma_{k=0}^{n}\left(\begin{array}{l} n \\ k \end{array}\right)y_{k} \leq \varepsilon . \end{aligned} $$

In MATLAB, this can be solved as follows

%Declaring input variables
n = 10; 
eps = 0.01;
delta = 0.2; 

for i = 1:n+1
    p(i, 1) = delta^(i-1)*(1-delta)^(n-i+1);
    nck(1,i) = nchoosek(n, i-1);
end
   
%Start optimization
cvx_begin        
   variable y(n+1)
   variable r(n+1)
   variable s
       
   minimize(s)
   subject to
        y >= 0;
        r >= 0;
        r <= s*ones(n+1, 1)
        y - r + p >= 0;
        nck*r ==1;
        nck*y <= eps;
cvx_end

How does one do vector constraints in Mathematica? Here is my attempt but I am not sure how to proceed from here

n = 10;
delta = 1/5;
eps = 1/50;
pvec = Array[p, n + 1];
nckvec = Array[nck, n + 1];
yvec = Array[y, n + 1];
rvec = Array[r, n + 1];

For[i = 1, i <= n + 1, i++, pvec[[i]] = delta^(i - 1)*(1 - delta)^(n - i + 1)]
For[i = 1, i <= n + 1, i++, nckvec[[i]] = Binomial[n, i]]

c1 = Thread[rvec . nckvec == 1];
c2 = Thread[yvec . nckvec <= eps];
c3 = Thread[yvec - rvec + pvec >= 0];
c4 = Thread[0 <= rvec];
c5 = Thread[0 <= yvec];
c6 = Thread[rvec <= s*ConstantArray[1, n + 1]];

constraints = Join[c1, c2, c3, c4, c5, c6];

I get an error when I try to join the constraints and I am also not sure how the final linear program solving using, say, LinearOptimization should be done.

$\endgroup$

1 Answer 1

6
$\begingroup$

Note: I have assumed that $p_k = \delta^k(1-\delta)^{n-k}$, as in OPs code.

Solution 1. One could use Minimize:

opt[delta_,eps_,n_]:=With[{p=(delta^#*(1-delta)^(n-#))&},
  Minimize[
    (*function and constraints*)
    {s,
    Table[y[k]-r[k]+p[k]>=0,{k,0,n}],
    Table[y[k]>=0,{k,0,n}],
    Table[0<=r[k]<=s,{k,0,n}],
    Sum[Binomial[n,k]*r[k],{k,0,n}]==1,
    Sum[Binomial[n,k]*y[k],{k,0,n}]<=eps},  
    (*variables*)
    Flatten[{s,Table[{y[k],r[k]},{k,0,n}]}]
]];

For example

opt[1/5,1/50,10]

gives

{1706527/19531250,{s->1706527/19531250,y[0]->0,r[0]->1706527/19531250,y[1]->1/500,r[1]->1126701/39062500,y[2]->0,r[2]->65536/9765625,y[3]->0,r[3]->16384/9765625,y[4]->0,r[4]->4096/9765625,y[5]->0,r[5]->1024/9765625,y[6]->0,r[6]->256/9765625,y[7]->0,r[7]->64/9765625,y[8]->0,r[8]->16/9765625,y[9]->0,r[9]->4/9765625,y[10]->0,r[10]->1/9765625}}

Solution 2. One could use LinearOptimization. It is essentially a copy and paste of solution 1 with some minor adjustments:

opt2[delta_,eps_,n_]:=With[{p=(delta^#*(1-delta)^(n-#))&},
  LinearOptimization[s,
    {Table[y[k]-r[k]+p[k]>=0,{k,0,n}],
     Table[y[k]>=0,{k,0,n}],
     Table[0<=r[k]<=s,{k,0,n}],
     Sum[Binomial[n,k]*r[k],{k,0,n}]==1,
     Sum[Binomial[n,k]*y[k],{k,0,n}]<=eps},
     Flatten[{s,Table[{y[k],r[k]},{k,0,n}]}]
]];

With this one should have more control, see for example the Method option.

Note. These solutions have the advantage that one can simply write down symbolic inequalities. Sometimes it may make sense to set up the problem in terms of matrices and use LinearProgramming.

$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.