x0 = [540 530];                            % initial guess
us = [30000 15000 15000 5000];             % operation condition
xs = trim('dens_io',x0',us',[],[],[1:4]);  % steady-state
[A,B,C,D]=linmod('dens_io',xs,us);         % linear model
x0 = xs;                                   % initial condition

