# LinearOptimization - how to solve with vector constraints?

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];
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.

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.