Here is a copy of my code:::

```
const = 40;
cvx_begin gp
variables b(n,1) epsilon
M = diag(b)*A + d;
Mhat = M;
for k = 1: r-1
Mhat = Mhat*M;
end
S2 = sum(Mhat,2);
maximize epsilon
subject to
S2 + epsilon*ones(n,1) <= const*ones(n,1);
sum(1/b2)<= 100;
b_lower*ones(n,1) <= b2 <= b_upper*ones(n,1);
cvx_end
```

It shows that my problem is solved, but the solution is infeasible. I checked it by plugging in the returned optimal b into M, and the first constraint is violated.