I am not sure how to treat the inputs tau and f into fmincon() because the objective function is not explicitly a function of tau or f.
Isn't it? In your formula for q_jk, all the tau's and f's appear on the right hand side. Anyway, the only thing that matters to fmincon is that you provide a function that will compute z for any given guess of the unknowns f and tau. It doesn't know or care what steps the function has to go through to reach z.
How do I model the inequality constraint that contains a summation?
That's what the Aineq,bineq arguments to fmincon are for.
Also, I need to satisfy this constraint for all i - can I do that in a single constraint?
No, you need n_i of them.
Can the solver handle two (one of which is multi-dimensional) variable solutions? The optimal f_ij array should be of size (i x j) and the optimal tau_j array should be of size (j x 1).
I would recommend, first of all, making a change of variables and replace exp(delta t/tau_j) with a single variable theta_j to simplify the computations and make them more linear. To answer your question though, f and tau (or f and theta) must be passed to your objective function bundled together in a single array unknowns. You can, of course, unpack them into separate variables within your objective function, if you wish.
To make it simpler to set this up, however, you can use some of the Optimization Toolbox tools for Problem-Based Optimization Setup together with the function prob2matrices() which must be downloaded. This will allow you to set up the optimization as below (where the blanks are for you to fill in):
f=optimvar('f',[ni,np],'LowerBound',0);
theta=optimvar('theta',np,'LowerBound',0,'UpperBound',1);
con.sumbound=sum(f,2)<=1;
sol0.f=____; %Initial guesses
sol0.theta=___;
[p,idx]=prob2matrices({f,theta},'Constraints',con,'x0',sol0);
fun=@(unknowns) zfun(unknowns,idx, I, qobj);
unknowns = fmincon(fun,p.x0,p.Aineq,p.bineq,[],[],p.lb,p.ub);
sol=x2sol(unknowns,idx);
f=sol.f; %final solutions
theta=sol.theta;
function z=zfun(unknowns,idx, I, qobj)
sol=x2sol(unknowns,idx);
f=sol.f;
theta=sol.theta;
Q=f.'*I;
q=nan(size(Q));
q(:,1)=_____;
for k=2:size(I,2)
q(:,k)=q(:,k-1).*theta + (1-theta).*Q(:,k);
end
z=norm(qobj-q,'fro').^2;
end
