In[]:=
(*TwoPeriod,TwoDepartmentModel*)​​c1=1;​​c2=2;​​r1=3;​​r2=6;​​delta=0.4;​​b=1;​​​​(*StationarityConstraint*)​​sc11=-c1+delta*r1-(-lambda11+mu1*(-r1)+v1*c1)==0;​​sc21=-c2+delta*r2-(-lambda21+mu1*(-r2)+v1*c2)==0;​​sc12=-delta^1*c1+delta^2*r1-(-lambda12+mu1*c1)==0;​​sc22=-delta^1*c2+delta^2*r2-(-lambda22+mu1*c2)==0;​​​​(*PrimaryFeasibility*)​​pf1={-a11≤0,-a12≤0,-a21≤0,-a22≤0};​​pf2=-(r1a11+r2a21)+(c1a12+c2a22)≤0;​​pf3=-b+(c1a11+c2a21)≤0;​​​​(*DualFeasibility*)​​df1={lambda11≥0,lambda12≥0,lambda21≥0,lambda22≥0};​​df2=mu1≥0;​​df3=v1≥0;​​​​(*ComplementarySlackness*)​​cs1={-lambda11*a11==0,-lambda12*a12==0,-lambda21*a21==​​0,-lambda22*a22==0};​​cs2=mu1*(-(r1*a11+r2*a21)+(c1*a12+c2*a22))==0;​​cs3=v1*(-b+(c1*a11+c2*a21))==0;​​​​(*OptimalPoints*)​​Reduce[Join[{sc11,sc21,sc12,sc22},pf1,{pf2,pf3},​​df1,{df2,df3},cs1,{cs2,cs3}],{a11,a12,a21,a22,lambda11,​​lambda12,lambda21,lambda22,mu1,v1},Reals]​​​​(*AnOptimalDiscountedTotalProfitValue*)​​b-(c1*a11+c2*a21)+delta*(r1*a11+r2*a21-(c1*a12+c2*a22))+​​delta^2(r1*a12+r2*a22)/.​​FindInstance[​​Join[{sc11,sc21,sc12,sc22},pf1,{pf2,pf3},df1,{df2,df3},​​cs1,{cs2,cs3}],{a11,a12,a21,a22,lambda11,lambda12,​​lambda21,lambda22,mu1,v1},Reals]
Reduce
:Reduce was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result.
Out[]=
0≤a11≤1.&&0≤a12≤3.&&a210.5(1.-1.a11)&&a220.5(3.a11-1.a12+6.a21)&&lambda110&&lambda120&&lambda210&&lambda220&&mu10.08&&v10.44
Out[]=
{1.44}
​
​
​
​
​
​