Initial code

.gitignore

@@ -0,0 +1,38 @@
+# Created by https://www.toptal.com/developers/gitignore/api/matlab
+# Edit at https://www.toptal.com/developers/gitignore?templates=matlab
+
+### MATLAB ###
+# Windows default autosave extension
+*.asv
+
+# OSX / *nix default autosave extension
+*.m~
+
+# Compiled MEX binaries (all platforms)
+*.mex*
+
+# Packaged app and toolbox files
+*.mlappinstall
+*.mltbx
+
+# Generated helpsearch folders
+helpsearch*/
+
+# Simulink code generation folders
+slprj/
+sccprj/
+
+# Matlab code generation folders
+codegen/
+
+# Simulink autosave extension
+*.autosave
+
+# Simulink cache files
+*.slxc
+
+# Octave session info
+octave-workspace
+
+# End of https://www.toptal.com/developers/gitignore/api/matlab
+

.session/default/.matlab/R2023a/MATLAB_Editor_State.xml

@@ -0,0 +1,24 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<Editor version="1.0">
+   <File absPath="/MATLAB Drive" lastWrittenTime="1680650292500" name="intersect.m">
+      <CodeFolds version="1.0"/>
+   </File>
+   <File absPath="/MATLAB Drive" lastWrittenTime="1680741646808" name="erlangtest.m">
+      <CodeFolds version="1.0"/>
+   </File>
+   <File absPath="/MATLAB Drive" lastWrittenTime="1680902189121" name="randomOreExtraction.m">
+      <CodeFolds version="1.0"/>
+   </File>
+   <File absPath="/MATLAB Drive" lastWrittenTime="1680903798551" name="solveforvolume.m">
+      <CodeFolds version="1.0"/>
+   </File>
+   <File absPath="/MATLAB Drive" lastWrittenTime="1680903938381" name="meanOreExtraction.m">
+      <CodeFolds version="1.0"/>
+   </File>
+   <File absPath="/MATLAB Drive" lastWrittenTime="1680981547995" name="findPriceQuantity.m">
+      <CodeFolds version="1.0"/>
+   </File>
+   <File absPath="/MATLAB Drive" lastWrittenTime="1681005657221" name="simulation.m">
+      <CodeFolds version="1.0"/>
+   </File>
+</Editor>

.session/default/.matlab/R2023a/VisibleSettings.json

@@ -0,0 +1 @@
+["Simulink","slhistory"]

.session/default/.matlab/R2023a/bookmarks.xml

@@ -0,0 +1 @@
+<?xml version="1.0" encoding="UTF-8"?><bookmarks><file xmlns="http://www.w3.org/1999/xhtml" absolutepath="/MATLAB Drive/simulation.m"></file><file xmlns="http://www.w3.org/1999/xhtml" absolutepath="/MATLAB Drive/intersect.m"></file><file xmlns="http://www.w3.org/1999/xhtml" absolutepath="/MATLAB Drive/findPriceQuantity.m"></file><file xmlns="http://www.w3.org/1999/xhtml" absolutepath="/MATLAB Drive/erlangtest.m"></file><file xmlns="http://www.w3.org/1999/xhtml" absolutepath="/MATLAB Drive/meanOreExtraction.m"></file><file xmlns="http://www.w3.org/1999/xhtml" absolutepath="/MATLAB Drive/randomOreExtraction.m"></file><file xmlns="http://www.w3.org/1999/xhtml" absolutepath="/MATLAB Drive/solveforvolume.m"></file></bookmarks>

.session/default/.matlab/R2023a/favorite_commands.json

@@ -0,0 +1 @@
+{"version":"1.1.1","state":"0","data":{"categories":{"CAT_1":{"label":"Favorite Commands","icon":"Category Icon","isInQAB":false,"showText":false,"editable":false,"state":"0"},"CAT_2":{"label":"Examples","icon":"Category Icon","isInQAB":false,"showText":false,"editable":true,"state":"0"}},"favorites":{"FAV_1":{"label":"About Favorite Commands","icon":"icon_help_favorite_16","code":"% Display help documentation for Favorite Commands\nhelpview('matlab', 'matlab_favorites');","isInQAB":false,"showText":false,"editable":true,"state":"0"},"FAV_2":{"label":"Clear Variables & Commands","icon":"fav_command_c","code":"% Remove all variables from the current MATLAB workspace\nclear;\n\n% Clear all input and output from the MATLAB Command Window\nclc;","isInQAB":false,"showText":false,"editable":true,"state":"0"},"FAV_3":{"label":"Go to User Folder","icon":"icon_favorite_command_16","code":"% Change MATLAB current folder to the first folder on the search path\ncd(userpath);","isInQAB":false,"showText":false,"editable":true,"state":"0"},"FAV_4":{"label":"MATLAB Logo","icon":"icon_matlab_favorite_16","code":"% Create a figure and display the MATLAB logo\nlogo;\n\n% Update figure so it is visible during animation\ndrawnow;\n\n% Rotate the MATLAB logo\n[az,el] = view;\nfor step = 1: 360\n    % Set new position of viewer\n    view(az + step, el);\n    % Stop MATLAB execution temporarily to slow down animation\n    pause(0.005);\nend\n","isInQAB":false,"showText":false,"editable":true,"state":"0"}}},"layout":{"categories":["CAT_1","CAT_2"],"favorites":{"CAT_1":[],"CAT_2":["FAV_1","FAV_2","FAV_3","FAV_4"]}}}

.session/default/.matlab/R2023a/matlab.prf

@@ -0,0 +1,10 @@
+#MATLAB Preferences
+#Sun Apr 09 02:04:54 GMT 2023
+RLHiddenB2_WB_2014b_1=I16376
+RLOrderB2_WB_2014b_1=S0:1:
+MatlabExitConfirm=Bfalse
+HtmlComponent.DefaultType.R2014b=SDUMMY
+CommandWindowClearConfirmation=Btrue
+RLWidthB2_WB_2014b_1_1=I130
+RLWidthB2_WB_2014b_1_0=I130
+RLPrevInitB2_WB_2014b_1=Btrue

.session/default/.matlab/R2023a/migratePref.txt

@@ -0,0 +1 @@
+T

.session/default/.matlab/R2023a/run_commands.m

@@ -0,0 +1,34 @@
+% MATLAB Editor Run Commands
+% @version 4
+% It is recommended that you do not edit this file directly, but use the
+% Run button. Changes made to comments that exist outside of a
+% configuration will not be saved.
+
+%% @name intersect
+%  @associatedFile /MATLAB Drive/intersect.m
+%  @mostRecentlyActioned true
+%  @uniqueId 7B2B0C0A
+
+intersect
+
+%% @name findPriceQuantity
+%  @associatedFile /MATLAB Drive/findPriceQuantity.m
+%  @mostRecentlyActioned true
+%  @uniqueId F5D8E1B1
+
+findPriceQuantity
+
+%% @name meanOreExtraction
+%  @associatedFile /MATLAB Drive/meanOreExtraction.m
+%  @mostRecentlyActioned true
+%  @uniqueId 7CF67C6B
+
+meanOreExtraction
+
+%% @name randomOreExtraction
+%  @associatedFile /MATLAB Drive/randomOreExtraction.m
+%  @mostRecentlyActioned true
+%  @uniqueId F426FBB3
+
+randomOreExtraction
+

.session/default/.matlab/R2023a/sl_toolstrip_plugins/preferences.json

@@ -0,0 +1 @@
+{"entries":[{"content":{"Toolstrip":{"content":{"FavoriteCommands":null,"QAB":null},"type":"prefs.Toolstrip","uuid":"b8d0f05f-b6f8-4a16-aef0-ce425512b1fc"}},"type":"prefs.Preferences","uuid":"14db1b6b-5aa2-46a7-b3c4-d972f73e8f75"}],"packageUris":["http://schema.mathworks.com/mf0/dig/1.0"],"version":"1.0"}

.session/default/.matlab/R2023a/sl_toolstrip_plugins/spritemap.dat

@@ -0,0 +1 @@
+d2fad784a7365cb60df1df475c5cbc7c

README.md

@@ -0,0 +1,132 @@
+# RimworldCommerceSimulator
+
+## Basic idea
+
+There are a series of villages that commerce between them.
+
+Assumptions:
+
+* There is a single good
+* There are consumer and producer villages
+* Each path between two villages adds a tranport cost to the good
+* The transport cost is proportional to the cost of the good
+* Each village extracts a random amount of gold each timestep
+* The offer price of the consumers increases with each gold extracted
+* Time is discrete (in days)
+
+## Simulator design
+
+### Village types
+
+Villages are divided into three types:
+
+* Producers: Produces the goods for a fixed cost each timestep
+* Consumers: Consumes a certain amount of goods each timestep
+* Traders: Try to learn the real equilibrium price for the good at each timestep, meanwhile it buys and sells goods
+
+### Villages and caravans
+
+Villages send out a caravan every T time, where T follows an Erlang distribution. The idea is that a caravan exiting the village is a [poisson point process](https://en.wikipedia.org/wiki/Poisson_point_process). The mean round trip time is given by the `distances` digraph:
+
+![Distances digraph](images/Transport.png)
+
+### Gold extraction
+
+Every village has raw gold waiting to be extracted. The gold extraction is done each timestep, with a Poisson distribution with mean given by the `meanOreExtraction` for each time t:
+
+![Gold mean](images/Gold_means.png)
+
+This function is $f(t) = a e^{bt} (1 - e^{-t})$. The parameters $b$ and $a$ are given by the system of equations:
+
+$$ max = \int_{0}^{\inf} f(t) dx = \frac{a}{b^2-b} $$
+
+$$ total = a e^{b ln(\frac{b-1}{b})} (1 - e^{-ln(\frac{b-1}{b})}) $$
+
+Where `max` is the maximum gold extraction in a day and `total` is the total extraction in a village's lifetime.
+
+Real accumulated extraction for each time step is given at random with Poisson distribution for each timestep, here is an example run:
+
+![Real accumulated extraction for each time step](images/Gold_total.png)
+
+For consumer villages, the gold extracted is used to increase the gold offered for each unit of the good. Also, the consumer villagers can never offer, more gold than is available at the village, if the offer is too high, it is updated maintaining the function slope.
+
+### Transport costs
+
+Traveling from village to village costs money, when a caravan travels from a village and back, it will try to buy goods from the other village, when buying it will use the `costs` digraph to calculate a "transport cost" that is proportional to the price payed for the goods:
+
+![Proportional costs](images/Costs.png)
+
+### Supply and demand and prices
+
+Each time a village sends a caravan, the offers and prices are equilized, and a new equilibrium price is calculated with the following equation:
+
+$$ \Delta price = \alpha (priceA - priceB) $$
+
+$$ priceA = priceA - \Delta price $$
+
+$$ priceB = priceB + \Delta price $$
+
+
+$$ \Delta offer = \alpha (offerA - offerB) $$
+
+$$ offerA = offerA - \Delta offer $$
+
+$$ offerB = offerB + \Delta offer $$
+
+Here is an example supply and demand graph:
+
+![Proportional costs](images/Price_and_offer_at_50_at_6.png)
+
+I call "offer" the amount the village thinks someone will pay for the good. I call "price" is the amount the village thinks someone will give up the good for. At the intersection is the guess for an equilibrium price.
+
+#### Trades
+
+A graph of the trades at each timepoint:
+
+![Trades](images/Trades.png)
+
+The average equilibrium price at each timepoint:
+
+![Goods at end](images/Prices.png)
+
+### Goods production
+
+At each timestep, the producer villages check if the offers are higher than the cost of producing the goods. It producers until the inventory is full or the offers are too low to cover the costs.
+
+## Results
+
+The simulation was run with 500 timesteps and the previously shown configuration. The producers are 1 and 8, the consumers are 3 and 4:
+
+### Aggregate gold at the end
+
+![Gold at end](images/Gold_per_village.png)
+
+### Aggregate goods at end
+
+![Goods at end](images/Goods_per_village.png)
+
+## Known issues
+
+### Traders will bankrupt themselves
+
+Trader villages have no memory, so will buy until their guessed equilibrium price says the price is too high. This means that when all the gold is extracted, they end up with lots of unsold goods that consumers have no gold for.
+
+Traders should probably use the newsvendor model:
+
+$$ q = F^{-1} (\frac{p-c}{p}) $$
+
+Although the goods are not perishible, the consumers are, this means that at some point the consumers with "run out", effectively destroying the goods value (see price graph).
+
+### Price update is somewhat broken
+
+The price equilibrium updating from village to village has problems keeping track of travel costs and the update function is too unrealistic.
+
+### The model does not support perishable goods
+
+This model is somewhat easily extended to perishable products. The only issue is that each perishable product has a memory, so they are not fungible anymore, this implies that when the product is near its expiration, the price decreases.
+
+The probability function for a good getting bad is: 
+
+$$ f(t) = min(S^{T-t},1) $$
+
+$T$ is the time when the product "dies". $S$ is the probability of selling the good in a single day.

findPriceQuantity.m

@@ -0,0 +1,29 @@
+function [equilibriumAmounts,equilibriumPrices] = findPriceQuantity(consumerOffers,producerPrices)
+%FINDPRICEQUANTITY Summary of this function goes here
+%   Detailed explanation goes here
+
+if max(producerPrices) < min(consumerOffers)
+    equilibriumPrices = (max(producerPrices)+min(consumerOffers))/2;
+    equilibriumAmounts = length(producerPrices);
+    return
+end
+
+aggreatePrices = reshape(producerPrices, 1, []);
+aggreatePrices = sort(aggreatePrices);
+
+aggreateOffers = reshape(consumerOffers, 1, []);
+aggreateOffers = sort(aggreateOffers, 'descend');
+
+equilibriumAmounts = [];
+equilibriumPrices = [];
+
+for i = 1:(min(length(aggreatePrices),length(aggreateOffers))-1)
+    [equilibriumAmount,equilibriumPrice] = intersect([i,i+1],aggreatePrices(i:i+1),[i,i+1],aggreateOffers(i:i+1));
+    equilibriumAmounts(i) = equilibriumAmount;
+    equilibriumPrices(i) = equilibriumPrice;
+end
+
+equilibriumAmounts(isnan(equilibriumAmounts)) = [];
+equilibriumPrices(isnan(equilibriumPrices)) = [];
+
+end

intersect.m

@@ -0,0 +1,37 @@
+function [outputX,outputY] = intersect(inputAX,inputAY,inputBX,inputBY)
+%INTERSECT Summary of this function goes here
+%   Detailed explanation goes here
+    x1 = inputAX(1);
+    x2 = inputAX(2);
+    x3 = inputBX(1);
+    x4 = inputBX(2);
+    y1 = inputAY(1);
+    y2 = inputAY(2);
+    y3 = inputBY(1);
+    y4 = inputBY(2);
+
+    A = [ det([x1 y1; x2 y2]) det([x1 1; x2 1])
+          det([x3 y3; x4 y4]) det([x3 1; x4 1])];
+    B = [ det([x1 1; x2 1]) det([y1 1; y2 1])
+          det([x3 1; x4 1]) det([y3 1; y4 1])];
+
+    outputX = det(A)/det(B);
+
+
+    C = [ det([x1 y1; x2 y2]) det([y1 1; y2 1])
+          det([x3 y3; x4 y4]) det([y3 1; y4 1])];
+    D = [ det([x1 1; x2 1]) det([y1 1; y2 1])
+          det([x3 1; x4 1]) det([y3 1; y4 1])];
+
+    outputY = det(C)/det(D);
+
+    if outputX < min([x1 x2 x3 x4]) || max([x1 x2 x3 x4]) < outputX
+        outputX = NaN;
+        outputY = NaN;
+    end
+
+    if outputY < min([y1 y2 y3 y4]) || max([y1 y2 y3 y4]) < outputY
+        outputX = NaN;
+        outputY = NaN;
+    end
+end

meanOreExtraction.m

@@ -0,0 +1,18 @@
+function [meanOre] = meanOreExtraction(time,maxExtraction,completeOreAvailability)
+%MEANOREEXTRACTION Summary of this function goes here
+%   Detailed explanation goes here
+
+% TODO solve analitically for a and b
+syms a b;
+M = maxExtraction;
+T = completeOreAvailability;
+
+S = vpasolve(a * exp(b * log((b-1)/b)) * (1-exp(-log((b-1)/b))) == M, a/(b*b-b) == T);
+clear a b;
+
+a = double(S.a);
+b = double(S.b);
+
+meanOre = a.*exp(b.*time).*(1-exp(-time));
+
+end

randomOreExtraction.m

@@ -0,0 +1,6 @@
+function [extractedOre] = randomOreExtraction(time,maxExtraction,completeOreAvailability)
+%RANDOMOREEXTRACTION Summary of this function goes here
+%   Detailed explanation goes here
+meanOre = meanOreExtraction(time,maxExtraction,completeOreAvailability);
+extractedOre = gamrnd(1,meanOre);
+end