Set S /1*2/; Parameters c(S) Driver cost /1 5,2 5/ A(s) Meeting function scaling parameter /1 3.25, 2 3.25/ alpha Meeting function elasticity /1 0.5, 2 0.5/ ; Scalar beta Customer value of time /27.69/ betab Background value of time /20/ kappa Elasticity /0.5/ kappab Background elasticity /0.01/ gamma Driver elasticity /0.25/ dr Delivery distance /4/ drb Background travel distance /5/ mu0 Consumer outside option /14.68/ omega0 Driver outside option /15/ ; Parameters l0 /2000/ N0 /800/ lb0 /20000/ v0 /30/ vc /0.0015/ ; Positive variables F(s) Fare for service s lambda(s) Demand for service s mu(s) Cost for service s wm(s) Meeting time for service s wr(s) Delivery time for service s dm(s) Meeting distance for service s v(s) Traffic speed N(s) Number of drivers for service s Ni(s) Number of idle drivers for service s omega(s) Hourly earning for service s R(s) Per ride earning for service s ; Variables PS(s) Firm s profit CS Consumer surplus BS Background surplus DS Driver surplus Z Objective function Nb ; alias(s2,s); Equations service_cost(s) service_demand(s) meet_time(s) meet_dist(s) service_time(s) service_fleet(s) back_fleet service_supply(s) fleet_cost(s) firm_profit(s) speed(s) obj consumer_surplus driver_surplus back_surplus ; *Bounds wm.up(s) = 1; wr.up(s) = 1; F.up(s) = 30; *Constraint definitions service_cost(s).. mu(s) =e= F(s) + beta*wm(s); service_demand(s).. lambda(s)*(sum(s2, exp(-kappa*mu(s2))) + exp(-kappa*mu0))=e= l0*exp(-kappa*mu(s)); meet_dist(s).. dm(s)*Ni(s)**alpha(s) =e= A(s); meet_time(s).. wm(s)*(v0 - vc*(sum(s2,N(s2))+Nb)) =e= dm(s); service_time(s).. wr(s)*(v0 - vc*(sum(s2,N(s2))+Nb)) =e= dr; service_supply(s).. N(s)*(sum(s2, exp(gamma*omega(s2))) + exp(gamma*omega0)) =e= N0*exp(gamma*omega(s)); service_fleet(s).. N(s) =e= Ni(s) + lambda(s)*(wm(s) + wr(s)); fleet_cost(s).. omega(s)*N(s) =e= R(s)*lambda(s) - c(s)*N(s); back_fleet.. Nb*(v0 - vc*(sum(s2,N(s2))+Nb))*exp(kappab*betab*drb/(v0 - vc*(sum(s2,N(s2))+Nb))) =e= lb0*drb; speed(s).. v0 - vc*(sum(s2,N(s2))+Nb) =e= v(s); firm_profit(s).. PS(s) =e= (F(s)-R(s))*lambda(s); consumer_surplus.. CS*kappa =e= log(sum(s2, exp(-kappa*mu(s2))) + exp(-kappa*mu0)); driver_surplus.. DS*gamma =e= log(sum(s2, exp(gamma*omega(s2))) + exp(gamma*omega0)); back_surplus.. BS*kappab*exp(kappab*betab*drb/(v0 - vc*(sum(s2,N(s2))+Nb))) =e= lb0*(1-(-kappab*betab*drb/(v0 - vc*(sum(s2,N(s2))+Nb)))) - kappab*betab*drb/(v0 - vc*(sum(s2,N(s2))+Nb))*lb0; obj.. Z =e= sum(s2,PS(s2)) + CS + DS + BS; model base /all/ duo /all-obj-back_surplus-driver_surplus-consumer_surplus/; *Solves base model to obtain starting values for JAMS option NLP = BARON; base.optfile = 1; Nb.lo = 0; solve base using NLP maximizing Z; *Uses EMP framework Nb.lo = -inf; file empinfo / '%emp.info%' /; put empinfo 'equilibrium' /; put 'implicit Nb back_fleet' /; loop(s, put 'max', PS(s), mu(s), lambda(s),F(s),R(s),omega(s),N(s),Ni(s),dm(s),wm(s),wr(s),v(s), service_cost(s), service_demand(s),meet_dist(s),meet_time(s),service_time(s),speed(s), service_supply(s),service_fleet(s),fleet_cost(s), firm_profit(s)/; ); putclose empinfo; *Option file for JAMS duo.optfile = 2; solve duo using emp;