Quelques calculs...

Apeiron Jeu 13 Juin - 11:16

Le but de ces deux pavés est de m'aider à réaliser le calcul demandé par Nimé pour implémenter le système économique et de lui permettre de comprendre et de vérifier la méthode.

Commençons par l'évolution de la population. Les règles nous disent :

Croissance démographique (en %) = 100*(1 - population actuelle/seuil démographique) + bonus des médicaments

Nouvelle population (en habitants) = (1 + croissance démographique/100)*ancienne population

A partir de maintenant, quand une de mes variables commencera par le caractère %, cela signifiera qu'elle est comprise entre 0 et 1. Je donne la présentation cette fois du point de vue du moteur du jeu et plus du point de vue utilisateur.

Nous avons donc :
croiss_pop = 1 - pop/seuil + %medic
pop = (1 + croiss_pop)*pop
d'où :
pop = (2 - pop/seuil + %medic)*pop

Cette formule a été conçue au départ pour être appliquée chaque jour (les jours sont donnés IRL) à partir d'une population initiale. Par exemple si le seuil est de 5, la population initiale de 1 et que la planète a des médicaments (disons +25% à la croissance) :
jour 0 : 1.000000 milliard
jour 1 : 2.050000 milliards
jour 2 : 3.772000 milliards
jour 3 : 5.641403 milliards
jour 4 : 6.328071 milliards
jour 5 : 6.229263 milliards
Après quelques oscillations ça se stabilise à exactement 6.25 milliards (5 + 25%) à partir du dixième jour.

Nous avons donc facilement les valeurs disons à minuit de la population. Le problème consiste à avoir la valeur en temps réel. Mais avant de le régler remarquons que même si on peut raisonnablement penser que seuil est une constante (qui ne sera modifiée que suite à une technologie spéciale ou à un changement de race du propriétaire), ce n'est pas du tout le cas de %medic, qui dépend de la richesse.

les médicaments donnent richesse*2,5% de bonus à la croissance démographique

Le +25% précédent correspondait donc à une richesse de 10. Notons que si la planète ne consommait pas de médicament alors la population et la richesse ne serait plus liées qu'indirectement pour la gestion des ressources.
Si la richesse est stabilisée bien sûr on peut considérer %medic comme constant, mais faisons le calcul en supposant que ce n'est pas encore le cas :

Croissance économique (en %) = (1 - richesse actuelle) + (100 - taux d'imposition)/10 + bonus des biens de consommation

Nouvelle population (en habitants) = (1 + croissance démographique/100)*ancienne population

croiss_rich = (1 - rich)/100 + (1 - %tax)/10 + %biens
rich = (1 + croiss_rich)*rich
d'où :
rich = (1 + (1 - rich)/100 + (1 - %tax)/10 + %biens)*rich

Il est à noter que cette formule ne dépend de la population que dans le sens où une évolution de la population déterminera une évolution du flux de production/consommation de biens. Je reviendrai sur les flux de ressource plus tard, pour l'instant supposons que les ressources sont infinies et calculons l'évolution de la richesse sur un exemple. Rappelons que :

les biens de consommation donnent richesse*0,5% de bonus à la croissance économique

Nous supposons ici que la richesse initiale est de 10 et que la taxation est constante à 20% :

jour 0 : 10.000000
jour 1 : 10.400000
jour 2 : 10.795199
jour 3 : 11.184086
jour 4 : 11.565235
jour 5 : 11.937333
jour 6 : 12.299193
jour 7 : 12.649770

Et cela nous donne l'évolution corrélée suivante pour la population :

jour 0 : 1.000000 milliard
jour 1 : 2.050000 milliards
jour 2 : 3.792500 milliards
jour 3 : 5.731909 milliards
jour 4 : 6.495516 milliards
jour 5 : 6.430741 milliards
jour 6 : 6.509744 milliards
jour 7 : 6.545750 milliards

La population ne continue d'évoluer que parce que la richesse et donc le bonus lié aux médicaments évolue. Contrairement à la population qui a tendance à se stabiliser au bout de quelques jours, la richesse évolue lentement.
A 20% de taxation la richesse converge vers 18, mais en partant de 10 elle n'est qu'à 12.65 au bout d'une semaine, 14.73 au bout de deux, 16.13 au bout de trois (c'est aussi là que la population passe à 7 milliards), 16.98 au bout de quatre, 17.46 au bout de cinq et 17.71 au bout de six (la population se stabilise vers 7.25 milliards, c'est à dire 5 + 18*2.5% comme attendu).
Evidemment, si on souhaite stabiliser la planète à 18 il vaut mieux partir d'une taxation de 0% et changer quand elle approchera la bonne valeur. Ainsi, à 0% en partant de 10, deux semaines suffisent pour atteindre 18, et un mois pour atteindre les 21 (limite de la richesse). Evidemment d'un point de vue gameplay on souhaiterait avoir une convergence plus rapide, mais d'un point de vue RP avoir une zone stable pendant longtemps pour atteindre le plafond économique ne semble pas déraisonnable. Cependant ce n'est pas le sujet.

Comme vu dans ce pavé, l'évolution de la population est corrélée à celle de la richesse, et les deux sont liées aux stocks. D'autres valeurs comme les revenus et les salaires sont également à considérer, mais dans le prochaine pavé je vais m'attacher à expliciter des méthodes d'approximation.

Apeiron Jeu 13 Juin - 12:19

La corrélation ayant été en partie traitée dans le précédent pavé, je vais maintenant m'attacher à répondre à la question : peut-on avec ces formules d'apparence simple, mais qui ne donnent que les valeurs à minuit, donner les valeurs intermédiaires ? Ceci étant motivé afin d'avoir un système économique en temps réel et non avec des tours d'une journée IRL.

J'ai commencé par évoquer la corrélation des différents paramètres afin de montrer que (si les ressources sont consommées [1]) il est nécessaire de calculer tous les paramètres avant de pouvoir donner les valeurs de minuits à minuits. Je vais supposer maintenant qu'elles sont connues et expliquer comment à partir de valeurs ponctuelles on peut tracer des valeurs intermédiaires raisonnables.

L'idée est de réussir à déterminer une fonction polynomiale passant par les points calculés : on appelle ça faire une interpolation polynomiale. Comme on a des valeurs très continues le résultat devrait être assez stable pour être raisonnable. Il faudra néanmoins discuter de l'erreur et de l'impact que cela peut avoir sur le jeu.

Les points sont calculés de minuit en minuit, aux dates t0, t1, t2, ... et on obtient comme dans le pavé précédent les valeurs de population p0, p1, p2, ... et de richesse r0, r1, r2, ... associées à nos pas de temps. Je vais faire l'exemple sur la population, la richesse se faisant exactement de la même manière. On cherche donc à obtenir un polynôme pop tel que pop(t0) = p0, pop(t1) = p1, pop(t2) = p2 ... à partir des points (t0, p0), (t1, p1), (t2, p2), ...

La méthode lagrangienne permet d'obtenir à partir de n+1 points un polynôme de degré inférieur ou égal à n. Par exemple si on ne prenait que l'état initial (t0, p0) et qu'on ne calculait que l'état du lendemain (t1, p1), on aurait un polynôme de degré 0 ou 1, c'est à dire une droite qui relierait les deux points. C'est l'approximation la plus stupide mais elle est séduisante car elle permettrait de gérer tous les calculs de façon linéaire. Discutons à présent de l'erreur.

Voici en gros la forme de la courbe qu'on devrait obtenir pour l'évolution de la population d'une planète fraichement colonisée : elle monte de plus en plus au début, puis la croissance se ralentit, la population dépasse un petit peu le maximum puis elle se stabilise. Je m'intéresse ici aux approximations entre t1 et t2, et entre t3 et t4 (en jaune) car les autres tronçons de courbe peuvent bien s'approximer par des droites : ce sont donc les pires cas.

Entre t1 et t2 l'approximation est bien au dessus de l'évolution continue : cela signifie que si un joueur mettait à jour les données de la planète (par exemple en changeant la taxation puis en annulant immédiatement le changement) il aurait des valeurs supérieures à la courbe rouge, et donc quand l'état serait recalculé à partir de la nouvelle valeur il aurait une nouvelle courbe rouge au dessus de la courbe actuelle. Réciproquement, entre t3 et t4 s'il mettait à jour il aurait des valeurs inférieures.

L'idée de l'interpolation polynomiale est qu'on n'atteindra jamais la courbe parfaite mais que plus on augmente le nombre de points dans une zone plus l'approximation devient précise. Ici on ne calcule pas plus de points entre t1 et t2, mais avec les points de t0 à t5 on obtient un polynôme de degré 5 passant par tous les points, en en raison de la régularité de nos valeurs l'évolution entre t1 et t2 devrait être plus conforme à une courbe comme la rouge. De plus, plus on plonge loin dans les calculs plus il devient difficile pour le joueur de theorycrafter à partir des valeurs intermédiaires obtenues. Evidemment, on peut aussi penser qu'un joueur capable de comprendre ce mécanisme peut mériter un gain en temps sur ces planète, ce gain étant de plus en plus réduit au fur et à mesure que la précision augmente.

Ainsi, il faudrait calculer loin dans l'avenir pour avoir des présents plus continus. Le problème c'est que le calcul en question a un coût. Avec n+1 points on obtient un polynôme de degré n, et pour le calculer il faut en gros inverser une matrice pleine (coût d'environ n^3). Il y a aussi des algorithmes plus optimisés comme ceux de Lagrange ou de Newton qui consistent à fabriquer des polynômes intermédiaires (la question est de savoir si le serveur manipule mieux les matrices ou les polynômes). Enfin, autant calculer plusieurs points dans l'avenir rendra la courbe plus propre (et donc moins arnaquable), autant je doute que calculer trop loin dans l'avenir soit bénéfique (typiquement la courbe de population a intérêt à s'arrêter à t5 ou t6 pour éviter de noyer les premières valeur par du bruit issu de la nécessité d'obtenir une courbe presque constante à partir de t6).

A ce stade il faut faire des essais numériques. Pour l'instant j'ai trouvé ce site qui utilise la méthode de Lagrange pour construire le polynôme. En le combinant avec ce site pouvant tracer des courbes j'obtiens pour notre exemple les courbes suivantes :

t0 : 1.000000 milliard
t1 : 2.050000 milliards
t2 : 3.792500 milliards
t3 : 5.731909 milliards
t4 : 6.495516 milliards
t5 : 6.430741 milliards
t6 : 6.509744 milliards
t7 : 6.545750 milliards

On part de l'état initial t0. La courbe bleue tient compte de tous les points jusqu'à t7, la rouge jusqu'à t6, la verte jusqu'à t5 et la grise jusqu'à t4 : on constate qu'essayer de trop tenir compte de l'avenir a introduit du bruit entre t0 et t1. Par contre les courbes sont remarquablement proches entre t1 et t3, et la grise décroche entre t3 et t4 (sa limite). Ainsi cet exemple numérique incite à chercher un polynôme de degré inférieur ou égal à (disons) 4. L'exemple sans médicament est similaire avec encore moins de divergences entre les courbes :

Remarquons que dans ces deux cas la forme générale de la courbe correspond bien aux prévisions, mais que les courbes sont très lisses, ce qui encourage davantage à prendre des polynômes de faible degrés, voire des droites. Intéressons-nous maintenant aux endroits courbés, typiquement entre t1 et t2 :

Et entre t3 et t4 :

La courbe est plus douce et probablement moins arnaquable en utilisant cette fois les approximations de plus haut degré.

Conclusion : Les degrés 1, 2 et 3 semblent trop mauvais en cas de courbes et les degrés 5, 6 et 7 semblent trop mauvais en cas de droite. Je pense donc utiliser l'interpolation polynomiale de degré 4. En partant sur ce postulat je ferai dans le prochain pavé des calculs explicites pour le serveur.

---------------------

[1] Un algorithme optimisé devra tenir compte des ordres de consommation afin de ne pas trop calculer quand les paramètres sont stables ou non corrélés.

Il est à noter que les stocks (de même que la taxation et le salaire) sont une intégrale quand la ressource est consommée par la population entière.

Apeiron Jeu 13 Juin - 16:56

Passons maintenant à la recherche de notre fameuse courbe rouge. Avant de passer à la partie explicite pour le serveur j'aimerais juste revenir sur les erreurs en montrant un enchainement de courbes de degré 4 :

La courbe bleue utilise les points t0 à t4.
La courbe rouge utilise les points t1 à t5.
La courbe verte utilise les points t2 à t6.
La courbe grise utilise les points t3 à t7.

Donc si on fait une mise à jour toutes les 24h, la bleue calcule entre t0 et t1, la rouge entre t1 et t2, la verte entre t2 et t3 et la grise entre t3 et t4. Nous constatons donc :
- la forme générale de la courbe est très bonne
- de légères discontinuités au moment des mises à jour
- les différences pas tout à fait négligeables n'apparaissent que quand on n'a pas mis à jour depuis plus de quelques jours
- comme les courbes se recroisent beaucoup (aux points utilisés pour les calculer, donc jusqu'à t + 4 jours) les erreurs sont automatiquement compensées sur des périodes de temps de quelques jours.

En somme, je ne vois que trois précautions à prendre :
- ne pas divulguer la méthode d'approximation : on approxime pour le confort de jeu, mais la méthode d'approximation ne fait pas partie du gameplay donc n'a aucune raison d'être publique, aussi prendre cette sécurité pourrait nous éviter des ennuis
- interdire les mises à jours trop rapprochées (difficile à estimer) afin d'éviter les discontinuités dont l'exploitation pourrait être faite par un bot avancé (le gars en question devrait cumuler mes compétences en maths et les compétences de Nimé en info), sans compter que cela pourrait gêner le serveur
- forcer des mises à jours de temps en temps (disons au moins tous les trois ou quatre jours) pour éviter de sortir de la fenêtre de sécurité.

A part cela, la méthode semble bonne.
De plus, grâce au faible degré du polynôme utilisé les calculs ne devraient pas être trop gourmands pour le serveur.

Je fatigue donc je posterai l'implémentation de l'algorithme plus tard.

Apeiron Ven 14 Juin - 8:38

Pour écrire l'algorithme il faut que je choisisse la méthode d'interpolation. En effet il est possible d'utiliser la méthode brute (consistant à inverser une matrice pleine ), lagrangienne (consistant à manipuler des polynômes formels) ou newtonienne (utilisant les différences divisées, méthode optimale pour un grand n). En effet, le n étant fixé petit il n'est pas sûr que la méthode newtonienne soit la meilleure.

Nous supposons déjà calculés les points (t0,y0), (t1,y1), (t2,y2), (t3,y3) et (t4,y4) car leur calcul devra être fait quelque soit la méthode d'interpolation choisie. L'idée est de calculer le polynôme pol tel que pour tout i pol(ti) = yi. pol(t) = a4*t^4 + a3*t^3 + a2*t^2 + a1*t + a0, donc pour calculer pol il faut calculer ses coefficients a4, a3, a2, a1 et a0.

Méthode brutale

Comme

a4*t0^4 + a3*t0^3 + a2*t0^2 + a1*t0 + a0 = y0
a4*t1^4 + a3*t1^3 + a2*t1^2 + a1*t1 + a0 = y1
a4*t2^4 + a3*t2^3 + a2*t2^2 + a1*t2 + a0 = y2
a4*t3^4 + a3*t3^3 + a2*t3^2 + a1*t3 + a0 = y3
a4*t4^4 + a3*t4^3 + a2*t4^2 + a1*t4 + a0 = y4

Cela revient à résoudre le système MA = Y, où M est la matrice :

1    t0    t0^2    t0^3    t0^4
1    t1    t1^2    t1^3    t1^4
1    t2    t2^2    t2^3    t2^4
1    t3    t3^2    t3^3    t3^4
1    t4    t4^2    t4^3    t4^4

où A est le vecteur :

a0
a1
a2
a3
a4

et Y le vecteur :

y0
y1
y2
y3
y4

Pour cela on applique le pivot de Gauss sur la matrice M|Y jusqu'à obtenir la matrice I|A. Le programme (optimisé pour tenir compte de la forme de la matrice, de tests inutiles et d'opérations inutiles à la fin sur M, mais il est illisible, navré : si tu veux l'algo regarde l'article Wikipedia) est le suivant :

for i = 1 to 4
   Y[i] = Y[i] - Y[0]
   for j = 0 to 4
      M[i,j] = M[i,j] - M[0,j]
for k = 1 to 4
   Y[k] = Y[k] / M[k,k]
   for j = k+1 to 4
      M[k,j] = M[k,j] / M[k,k]
   M[k,k] = 1
   for i = k+1 to 4
      Y[i] = Y[i] - M[i,k]*Y[k]
      for j = k+1 to 4
         M[i,j] = M[i,j] - M[i,k]*M[k,j]
      M[i,k] = 0
(à ce stade la matrice a des 1 sur la diagonale et des 0 en dessous, il ne reste qu'à remonter le système : )
for k = 3 to 0 (seule boucle du programme à décrémenter)
   for i = k+1 to 4
      Y[k] = Y[k] - M[k,i]*Y[i]
return Y

car on a maintenant :

a0 = Y[0]
a1 = Y[1]
a2 = Y[2]
a3 = Y[3]
a4 = Y[4]

Ensuite, pour calculer pol(t) il suffit de faire pol(t) = a4*t^4 + a3*t^3 + a2*t^2 + a1*t + a0, ce qui rajoute 1 = , 10 * et 4 +

Le coût total est de :
245 =
4 +
202 -
188 *
28 /
(si je ne me suis pas planté, parce que ça parait beaucoup... si je me suis planté c’est probablement d'un facteur 6)

Méthode de Lagrange

Les deux méthodes suivantes consistent en gros à gratter de façon ingénieuse pour avoir une matrice qui s'inverse mieux, dans les deux cas en représentant pol(t) dans une autre base de polynômes que 1, t, t^2, t^3 et t^4.

Ici : pol(t) = y0*L0(t) + y1*L1(t) + y2*L2(t) + y3*L3(t) + y4*L4(t), où :
L0(t) = (t - t1)(t0 - t1) + (t - t2)(t0 - t2) + (t - t3)(t0 - t3) + (t - t4)(t0 - t4)
L1(t) = (t - t0)(t1 - t0) + (t - t2)(t1 - t2) + (t - t3)(t1 - t3) + (t - t4)(t1 - t4)
L2(t) = (t - t0)(t2 - t0) + (t - t1)(t2 - t1) + (t - t3)(t2 - t3) + (t - t4)(t2 - t4)
L3(t) = (t - t0)(t3 - t0) + (t - t1)(t3 - t1) + (t - t2)(t3 - t2) + (t - t4)(t3 - t4)
L4(t) = (t - t0)(t4 - t0) + (t - t1)(t4 - t1) + (t - t2)(t4 - t2) + (t - t3)(t4 - t3)

On obtient une expression directe mais compliquée pour pol(t), contrairement à la méthode brutale où il fallait faire pas mal de calcul en amont mais où on obtenait directement les coefficients de pol. Comme n est petit, la méthode brutale aurait pu être meilleure si on appelait plusieurs fois la fonction pol obtenue, mais en réalité on ne la calcule que pour l'utiliser une fois avant de mettre à jour, donc la méthode lagrangienne est ici bien supérieure.

Plus précisément le coût est de :
5 =
19 +
40 -
25 *

(on peut éventuellement précalculer certains termes pour échanger des - contre des =)

Méthode de Newton

Enfin cette méthode consiste à trouver une meilleure base de polynôme en utilisant les différences divisées.
Ici : pol(t) = y0 + [y0,y1](t - t0) + [y0,y1,y2](t - t0)(t - t1) + [y0,y1,y2,y3](t - t0)(t - t1)(t - t2) + [y0,y1,y2,y3,y4](t - t0)(t - t1)(t - t2)(t - t3), où les [y...] sont les différences divisées calculées ainsi :
[y0,y1] = (y1 - y0)/(t1 - t0)
[y1,y2] = (y2 - y1)/(t2 - t1)
[y2,y3] = (y3 - y2)/(t3 - t2)
[y3,y4] = (y4 - y3)/(t4 - t3)
[y0,y1,y2] = ([y1,y2] - [y0,y1])/(t2 - t0)
[y1,y2,y3] = ([y2,y3] - [y1,y2])/(t3 - t1)
[y2,y3,y4] = ([y3,y4] - [y2,y3])/(t4 - t2)
[y0,y1,y2,y3] = ([y1,y2,y3] - [y0,y1,y2])/(t3 - t0)
[y1,y2,y3,y4] = ([y2,y3,y4] - [y1,y2,y3])/(t4 - t1)
[y0,y1,y2,y3,y4] = ([y1,y2,y3,y4] - [y0,y1,y2,y3])/(t4 - t0)

Le coût est de :
11 =
4 +
30 -
10 *
10 /

Conclusion

Cela se joue entre Lagrange :
5 =
19 +
40 -
25 *

Et Newton :
11 =
4 +
30 -
10 *
10 /

Si jeune m'abuse les = sont négligeables et les + comptent comme des -, donc à moins que les / soient plus couteux que des * il vaut mieux Newton. Sinon, Lagrange ira très bien.

C'est à toi de trancher, Nimé ! Smile

Je te rédigerai l'algo complet après.

Nimeroni Ven 14 Juin - 15:07

Bien que le coût exact dépend du processeur, la division est généralement très coûteuse.

Voici un exemple pour un processeur grand publique, le Pentium, avec des valeurs entières:
=, +, - : 1 cycle
* : 9 cycles
/ : de 17 / 25 / 41 cycles (pour 8 / 16 / 32 bits)

Note: si tu multiplie ou que tu divise par une puissance de 2, le coût passe a 1 cycle.

Apeiron Ven 14 Juin - 17:15

Donc ça fait pour Lagrange :
5 + 19 + 40 + 9*25 = 289

Et pour Newton :
11 + 4 + 30 + 9*10 + 17*10 = 304 (au mieux, mais comme on manipule de grands entiers...)

Va pour du Lagrange. Smile

PS : En fait non, je me suis planté en écrivant les polynômes de Lagrange. Je referai les calculs et les vérifications plus tard.

Apeiron Ven 14 Juin - 18:16

En fait pour Lagrange : pol(t) = y0*L0(t) + y1*L1(t) + y2*L2(t) + y3*L3(t) + y4*L4(t), où :
L0(t) = (t - t1)/(t0 - t1) * (t - t2)/(t0 - t2) * (t - t3)/(t0 - t3) * (t - t4)/(t0 - t4)
L1(t) = (t - t0)/(t1 - t0) * (t - t2)/(t1 - t2) * (t - t3)/(t1 - t3) * (t - t4)/(t1 - t4)
L2(t) = (t - t0)/(t2 - t0) * (t - t1)/(t2 - t1) * (t - t3)/(t2 - t3) * (t - t4)/(t2 - t4)
L3(t) = (t - t0)/(t3 - t0) * (t - t1)/(t3 - t1) * (t - t2)/(t3 - t2) * (t - t4)/(t3 - t4)
L4(t) = (t - t0)/(t4 - t0) * (t - t1)/(t4 - t1) * (t - t2)/(t4 - t2) * (t - t3)/(t4 - t3)

Ainsi, on a bien Li(tj) = 1 si i = j et 0 sinon, d'où pol(ti) = yi.
On a donc un coût de :
6 =
4 +
40 -
15 *
20 /

Contre :
11 =
4 +
30 -
10 *
10 /

Conclusion : Newton est bien meilleur.

Apeiron Ven 14 Juin - 19:12

// Variables utilisées :
// float pop (en milliards)
// int seuil (en milliards)
// float richesse (entre 1 et 21 en gros)
// float tax (entre 0 et 1)
// t (date actuelle)

// Base de données :
float medic = 0.025
float biens = 0.005

// Calcule l'état de la variable choisie en t à partir des points (t0,y0), (t1,y1), (t2,y2), (t3,y3) et (t4,y4) déjà calculés
pol(t, t0, t1, t2, t3, t4, y0, y1, y2, y3, y4) =
   y01 = (y0 - y1)/(t0 - t1)
   y12 = (y1 - y2)/(t1 - t2)
   y23 = (y2 - y3)/(t2 - t3)
   y34 = (y3 - y4)/(t3 - t4)
   y012 = (y01 - y12)/(t0 - t2)
   y123 = (y12 - y23)/(t1 - t3)
   y234 = (y23 - y34)/(t2 - t4)
   y0123 = (y012 - y123)/(t0 - t3)
   y1234 = (y123 - y234)/(t1 - t4)
   y01234 = (y0123 - y1234)/(t0 - t4)
   r = y0 + y01(t - t0) + y012(t - t0)(t - t1) + y0123(t - t0)(t - t1)(t - t2) + y01234(t - t0)(t - t1)(t - t2)(t - t3)
   return r

// On fait la mise à jour en t en souhaitant connaître l'état actuel de la planète en sachant qu'on connait déjà
// les valeurs prévues lors de la dernière mise à jour pour la richesse et la population, stockées dans les tableaux R et T
// initialisation : richesse et pop valent 1 si colonisation, valeurs de l'ancien propriétaire sinon
update(t,t0,R,T) =
   richesse = pol(t, t0, t0+1j, t0+2j, t0+3j, t0+4j, R[0], R[1], R[2], R[3], R[4])
   pop = pol(t, t0, t0+1j, t0+2j, t0+3j, t0+4j, P[0], P[1], P[2], P[3], P[4])
   R[0] = richesse
   P[0] = pop
   for i = 0 to 3
      R[i+1] = (1 + (1 - R[i])/100 + (1 - tax)/10 + R[i]*biens)*R[i]
      P[i+1] = (2 - P[i]/seuil + R[i]*medic)*P[i]
   t0 = t
   t = newupdate(planète)
   return (t,t0,R,T)

// newupdate dépend passivement non seulement de R et P, mais aussi de tous les déclencheurs liés aux stocks,
// aux crédits, à la recherche, aux actions des joueurs, etc.

// On a ainsi pour chaque planète une transition d'état (t,t0,R,T) = update(t,t0,R,T)

Apeiron Ven 14 Juin - 20:22

Bon, à moins que tu aies une demande particulière, je ne vais pas m'avancer davantage pour l'algorithme de mise à jour. Après tout, c'est ton domaine. Par contre, maintenant que j'ai réglé le calcul de la population et la richesse courantes je vais pouvoir m'attaquer aux calculs secondaires, par exemple :

Crime (en %) = distance (4 max) *5 + 20-richesse
taxation = taux d'imposition*richesse de la planète*population actuelle (en milliards)
salaires = pop_active*richesse (sans compter le trait financier)

On ne peut pas se contenter de les calculer en utilisant la valeur actuelle de population et de richesse : ils vont évoluer en même temps qu'elles, et donc il faut évaluer l'effet cumulé dans le temps. Il nous faut donc calculer des intégrales :

taxation[t0 to t1] = tax*integ[t = t0 to t1](richesse(t)pop(t)dt)
salaires[t0 to t1] = pop_active*integ[t = t0 to t1](richesse(t)dt)

comme richesse(t) et pop(t) sont approximées par des polynômes il est possible de calculer leur primitive, mais pour cela il faut extraire explicitement les valeurs des coefficients, ce qui va changer un peu l'algo mais qui devrait rester faisable vu que n = 4.

A ce stade, ce qui m'inquiète davantage c'est de pouvoir détecter un changement d'état (notamment baisse de population active, rupture de stock pour une ressource consommée par la population et banqueroute financière). En effet, ce n'est pas la même chose de demander la valeur à une date et de deviner la date où la valeur s'annule.

Apeiron Sam 15 Juin - 16:40

C'est délicat comme problème : les méthodes d'approximations requièrent généralement des propriétés sur la courbe pour assurer un résultat. Du coup je réfléchis à plusieurs pistes :
- comprendre pourquoi la richesse est en pratique une droite [1], et en déduire d'éventuelles simplifications
- esquisser le fait qu'à part les variations quand la planète finit par se stabiliser les quantités regardées sont strictement monotones [2]
- voir si du coup des algos pas trop couteux peuvent être trouvés, et sinon redéfinir le système économique sous une forme plus propice aux mises à jour en tenant compte de ce que cette analyse m'a appris.

[1] Bon, presque (courbe en bleu, droite en rouge) :

En fait c'est juste l'évolution sur quatre jours, les ordonnées étant entre 10 et 11. Il est probable que ça ressemble aussi violemment à une droite parce que l'évolution est très lente, contrairement à la population.
A titre de comparaison, avec toutes les ressources disponibles une planète de seuil 7 nouvellement colonisée atteint le seuil planétaire en 5 jours (et après suit l'évolution de la richesse), alors qu'il faut 27 jours pour atteindre 10 de richesse et 56 pour atteindre la maximum de 21. Autant il peut être intéressant que des valeurs prennent du temps à se mettre en place pour donner de l'importance aux zones pacifiées, autant là c'est long je trouve.
PS : Augmenter la valeur de la constante liée aux biens augmente le plafond mais pas la convergence, et la monter à 1% et plus entraine une divergence.

[2] Par exemple pour la population :

PS : voici l'évolution de la richesse sur 50 jours (j'ai affiché la droite pour mieux montrer la courbure) :

Nimeroni Sam 15 Juin - 21:33

2 mois, c'est un peu trop. Qu'est-ce qui augmente la convergence alors ?

Apeiron Sam 15 Juin - 23:31

Je planche sur un système alternatif davantage basé sur les mises à jours, utilisant la fonction log comme oracle.

On souhaite qu'à partir d'une date initiale t0 on ait le même résultat en t2 que si entretemps on fait une mise à jour en t1.
Supposons qu'on calcule la population selon la formule : pop(t2) = pop(t0) + log(t2/t0).
Du coup, quand on fait la mise à jour à un état intermédiaire, on a :
pop(t1) = pop(t0) + log(t1/t0) et pop(t2) = pop(t1) + log(t2/t1), d'où :
pop(t2) = pop(t1) + log(t2/t1) = pop(t0) + log(t1/t0) + log(t2/t1) = pop(t0) + log(t1/t0 * t2/t1) = pop(t0) + log(t2/t0)

Ainsi, si l'accroissement de population entre t0 et t1 est donné par log(t1/t0) on a une formule explicite et stable par mise à jour. Pourquoi le log et pas l'exponentielle ou une droite ? Parce qu'il correspond à l'idée d'avoir une forte croissance au début et une faible à la fin, ce qui permet de rapidement faire des choses avec ses planètes, tout en devant attendre pour être à un niveau optimal :

Quelques calculs... Logarithms

(ne regarder que la partie positive des courbes bien sûr)
Évidemment, comme la fonction log tend vers l'infini, il faudrait poser que la population ne peut dépasser le seuil de la planète, ce qui est en fait plus habituel que l'ancien système.

Comment gérer l'évolution de la croissance, et notamment la décroissance ?
log est le logarithme binaire (j'ai supposé qu'il serait plus naturel pour le serveur), ln est le logarithme népérien (fonction rouge), log10 est le logarithme en base décimale (fonction verte). et plus généralement loga est le logarithme en base a On passe d'une base à l'autre par la formule :
logb(x) = loga(x)/loga(b)
Réciproquement, cela signifie que si c est une constante différente de 0 alors log(x)/c est encore une fonction logarithmique donc stable par mise à jour. Si c est positif on a une croissance, et sinon on a une décroissance. Ce c deviendra donc notre modificateur. Comment le choisir ?
log(2) = 1, ln(e) = 1, log10(10) = 1, et plus généralement loga(a) = 1. Supposons pour l'instant que 1 est le plafond, cela signifie que si on veut l'atteindre à la date tmax il faut utiliser loga(x) = log(x)/log(a). Donc, en partant d'une population initiale nulle [1] il faudrait un modificateur égal à log(tmax) pour atteindre 1 en tmax. Évidemment, si on veut un plafond de popmax au lieu de 1, il suffit d'utiliser popmax*log(x)/log(a).

[1] log suppose implicitement un t0 fixe pour une population initiale donnée. Par exemple on peut partir d'une population nulle avec t0 = 1, mais cela signifiera que les mises à jours se font non pas sur la date en cours mais plutôt en un sens sur l'âge de la planète. Je développerai ceci plus clairement avec les exemples numériques.

Apeiron Dim 16 Juin - 0:50

ça devrait donner quelque chose dans ce goût-là :

// base de données et valeurs dérivées
popInit = 0
ageInit = exp(popInitial)
bonusMedic = 0.025
tMax = 1 mois (en secondes)
tauxPop = log(tMax)

// initialisation des valeurs à la colonisation / conquête
colo(joueur, planete, t0) =
   planete.proprio = joueur
   planete.seuil = seuilPop(joueur, planete)
   planete.pop = popInit
   planete.age = ageInit
   planete.upd = t0

// mise à jour de la population en t1 en supposant t0 connu
update(planete, t1) =
   joueur = planete.proprio
   plafond = planete.seuil
   if(joueur.tech(biosphere)):
      plafond = plafond + 1
   if(planete.conso(medic)):
      plafond = plafond * (richesse * bonusMedic)
   population = planete.pop
   if population != plafond:
      t0 = planete.upd
      a0 = planete.age
      a1 = a0 + (t1 - t0)
      deltaPop = plafond*(log(a1/a0)/tauxPop)
      if population < plafond:
         if deltaPop > (plafond - population):
            deltaPop = plafond - population
         population = population + deltaPop
      if population > plafond:
         if deltaPop > (population - plafond):
            deltaPop = population - plafond
         population = population - deltaPop
      planete.pop = population
      planete.age = exp(population)
   planete.upd = t1

Apeiron Dim 16 Juin - 13:00

Code:: int main() { float pop = 0.0; float seuil = 7.5; float taux = log(30); int jour = 1; int pas = 1; float a0, a1; char c; while(1){ printf("jour : %d\n",jour-1); printf("pop : %f\n",pop); a0 = jour; a1 = jour + pas; pop = pop + seuil*(log(a1/a0)/taux); if(pop > seuil){pop = seuil;} scanf("%c", &c); jour = a1; } return 0; }

Les résultats sont bien invariants selon la taille des pas que je prends.
On a bien une croissance rapide au début et faible à la fin.

Plus précisément, pour une pop initiale de 0, un seuil de 5 et +50% grâce aux medocs, il faut :
- 1 jour pour atteindre 1 milliard
- 2 jours pour atteindre 2 milliards
- 3 jours pour atteindre 3 milliards
- 6 jours pour atteindre 4 milliards
- 9 jours pour atteindre 5 milliards
- 15 jours pour atteindre 6 milliards
- 23 jours pour atteindre 7 milliards
- 30 jours pour plafonner à 7.5 milliards.

La croissance démographique fonctionne donc correctement. Smile

Faut que je teste maintenant la décroissance et des mises à jours sur les états en cours de route.

Apeiron Dim 16 Juin - 18:11

Code:: int main() { float pop = 9.0; float seuil = 0.0; float taux = log(30.0); float age = 1.0; float pas = 1.0; float a0, a1, delta; char c; while(1){ printf("t : %f\n",age); printf("pop : %f\n",pop); a0 = age; a1 = age + pas; if(seuil > pop){ delta = seuil*(log(a1/a0)/taux); if(delta > (seuil - pop)){ delta = seuil - pop; } pop = pop + delta; } if(pop > seuil){ delta = 9.0*(log(a1/a0)/taux); if(delta > (pop - seuil)){ delta = pop - seuil; } pop = pop - delta; } scanf("%c", &c); age = a1; } return 0; }

Supposons que pour une raison ou une autre (exemple : pandémie Very Happy

) la population décroisse jusqu'à extinction :

(population initiale de 9 milliards d'habitants, extinction en 30 jours)

La partie décroissante du programme contient la valeur de la population initiale, il faut que je règle ceci dans le cas général, puis que j'implémente une version basée non pas sur l'âge mais sur la date.

Apeiron Dim 16 Juin - 23:32

Je pense être arrivé à gérer les cas continus de croissance et de décroissance, avec des valeurs de population initiale et de plafond quelconques.

Code:: int main() { float pop0 = 6.0; float seuil = 9.0; float pop = pop0; float time = 30.0; float taux = log(time); float pas = 1.0; float a0, a1, delta; char c; if(seuil > pop){ a0 = exp(taux*pop0/seuil); } if(seuil < pop){ a0 = 1; } while(1){ printf("t : %f\n",a0); printf("pop : %f\n",pop); a1 = a0 + pas; if(seuil > pop){ delta = seuil*(log(a1/a0)/taux); if(delta > (seuil - pop)){ delta = seuil - pop; } pop = pop + delta; } if(pop > seuil){ delta = pop0*(log(a1/a0)/taux); if(delta > (pop - seuil)){ delta = pop - seuil; } pop = pop - delta; } scanf("%c", &c); a0 = a1; } return 0; }

Exemple 1

Cas d'une croissance à partir d'une population initiale :
pop0 = 6
seuil = 9

(oui, après 10 jours la courbe se confond avec la croissance à partir d'une population initiale de 0)

Exemple 2

Cas d'une décroissance jusqu'à un nouveau seuil :
pop0 = 9
seuil = 3

(oui, la courbe se confond avec la décroissance jusqu'à 0 pour la même population initiale, et ce jusqu'à atteindre la valeur de seuil)

Bien sûr j'ai vérifié à chaque fois une invariance par dilatation du pas de temps.

La prochaine étape sera de tracer des courbes avec des changements de paramètres à la volée.

Apeiron Lun 17 Juin - 12:19

La variable age sert à tenir compte du fait que l'état actuel n'est probablement pas un état qui était stabilisé.

Par exemple si on part d'une population initiale de 0 et qu'on désire atteindre la valeur de seuil en 30 jours, on aura que la population à l'instant 1 < t < 30 est égale à :
pop(t) = seuil*log(t)/log(30)
(on a bien pop(1) = 0 et pop(30) = seuil)

A contrario, si on part d'une population initiale de pop0 il faudra moins de temps pour atteindre le même seuil s0, et la croissance doit être la même que sur la courbe précédente à l'instant t où la population atteint pop0. Donc il ne faut pas commencer à 1 mais à un âge a0 vérifiant :
pop0 = s0*log(a0)/log(30)
d'où a0 = exp(log(30)*pop0/s0).

Le problème est le même que pour la décroissance. Si on part d'une population stabilisée à l'ancien seuil et qu'on essaie de la réduire à 0 en 30 jours alors la population à l'instant 0 < t < 30 est égale à :
pop(t) = seuil - seuil*log(t)/log(30)
(on a bien pop(1) = s0 et pop(30) = 0)

Donc, en partant d'une population initiale pop0 et d'un seuil initial supérieur s0, la planète doit avoir un âge a0 vérifiant :
pop0 = s0 - s0*log(a0)/log(30)
d'où a0 = exp(log(30)*(1 - pop0/s0))

Cela règle mon problème lors d'une décroissance à population initiale non stabilisée : la décroissance prenait la population actuelle comme valeur de référence donc faire une mise à mise aurait accru la décroissance, ce qui aurait rompu l'invariance. Il fallait mémoriser la population stabilisée précédente, mais celle-ci aurait très bien pu n'avoir jamais été atteinte. C'est pour cela que la valeur à mémoriser n'était pas l'ancienne population mais l'ancien seuil.

Du coup, quand une décroissance a commencé (par exemple parce que la planète n'a plus de médicament), le joueur aurait pu auparavant redonner des médicaments puis à l'instant d'après les retirer, ce qui aurait mis à jour l'ancienne population par la population actuelle. Maintenant on avait auparavant un ancien seuil avec médicaments quand la décroissance commence, et quand le joueur donne puis reprend les médicaments on a de nouveau un ancien seuil avec médicaments, donc la décroissance continu au même rythme.

Conclusion : quand on met à jour l'état de la planète, la nouvelle valeur de population est obtenu à partir de l'ancienne, de l'ancien seuil et du nouveau seuil. Le calcul est fait en utilisant l'état d'avancement de la planète, qu'on calcule à partir de ces données et de la durée depuis la dernière mise à jour.

Apeiron Mar 18 Juin - 14:48

Maintenant l'affichage se fait selon des dates réelles et non l'âge, et je pense avoir réglé mon problème quant à la décroissance relative suite aux calculs du dernier post. Il ne faut pas mémoriser l'âge mais le dernier seuil calculé, ce qui est fait de façon brutale dans le programme suivant :

Code:: int main(){ float pop = 0.0; float s0, s1; float t0 = 10.0; float tmax = 30.0; float pas = 1.0 ; float time = 30.0; float taux = log(time); float t1, a0, a1, s, delta; printf("date : %f\n",t0); printf("pop : %f\n\n",pop); for (t1 = t0 + pas; t1 <= tmax; t1 = t1 + pas){ if(t1 <= 27){ s0 = 0.0; s1 = 9.0; } else{ s0 = 9.0; s1 = 6.0; } if(s1 > pop){ a0 = exp(taux*pop/s1); a1 = a0 + (t1 - t0); delta = s1*(log(a1/a0)/taux); if(delta > (s1 - pop)){delta = s1 - pop;} pop = pop + delta; } if(pop > s1){ a0 = exp(taux*(1 - pop/s0)); a1 = a0 + (t1 - t0); delta = s0*(log(a1/a0)/taux); if(delta > (pop - s1)){delta = pop - s1;} pop = pop - delta; } t0 = t1; printf("date : %f\n",t0); printf("pop : %f\n\n",pop); } return 0; }

L'exemple suivant montre une planète de seuil 6 colonisée au jour 10, avec des médicaments dont l'approvisionnement est suspendu au jour 27. Le joueur fait ici une mise à jour tous les jours :

(j'ai bien sûr vérifié l'invariance par mise à jour ainsi que les croissances ou décroissances asymptotiques)

L'exemple suivant correspond à une population de 9 milliards soudain touchée par la maladie (le seuil est ajusté à 0) au jour 10, et le remède est trouvé au tour 20 :

La forme de la courbe est bonne, mais j'ai découvert que le seuil de population est de nouveau atteint le jour 50, ce qui signifie que le taux de croissance est calculé selon la population de survivants et non selon une population de 0 comme cela devrait l'être... je ne comprends pas pourquoi...

Edit :

Code:: int main(){ // population initiale, les seuils sont fixés dans le for float pop = 9.0; float s0, s1; // l'analyse est faite de t0 à tmax float t0 = 10.0; float tmax = 50.0; float pas = 1.0 ; // time est la durée pour passer du min au max float time = 30.0; float taux = log(time); float t1, a0, a1, s, delta; printf("date : %f\n",t0); printf("pop : %f\n\n",pop); for (t1 = t0 + pas; t1 <= tmax; t1 = t1 + pas){ // mise à jour des seuils à certaines dates if(t1 <= 20){ s0 = 9.0; s1 = 0.0; } else{ s0 = 0.0; s1 = 9.0; } // croissance de la population if(s1 > pop){ a0 = exp(taux*pop/s1); a1 = a0 + (t1 - t0); delta = s1*(log(a1/a0)/taux); if(delta > (s1 - pop)){delta = s1 - pop;} pop = pop + delta; } // décroissance de la population if(pop > s1){ a0 = exp(taux*(1 - pop/s0)); a1 = a0 + (t1 - t0); delta = s0*(log(a1/a0)/taux); if(delta > (pop - s1)){delta = pop - s1;} pop = pop - delta; } // mise à jour de l'état t0 = t1; printf("date : %f\n",t0); printf("pop : %f\n\n",pop); } return 0; }

En fait il ne s'arrête pas à 50 mais à 48, or 2 jours c'est bien le temps nécessaire pour atteindre la population de survivants si on était parti d'une population de 0. Donc en fait tout va bien. Very Happy

Apeiron Mer 19 Juin - 14:13

L'algo pour l'évolution de la population est donc :

Code:: data structure = time = 2592000 // 30 jours en secondes rate = log(time) bonusMedic = 0.5 // +50% au seuil struct planet = size climate owner pop thresh lastThresh lastUpdate popThresh(planet) = thresh = planet.size if not habitable(planet.climate, planet.owner) : thresh = 0 if not favorable(planet.climate, planet.owner) : thresh = thresh / 2 if conso(planet, medic) : thresh = thresh * (1 + bonusMedic) return thresh colo(planet, player, t) = planet.owner = player planet.pop = 0 planet.thresh = popThresh(planet) planet.lastThresh = 0 planet.lastUpdate = t popUpdate(planet,t1) = s0 = planet.lastThresh s1 = planet.thresh s = popThresh(planet) if s != s1 : s0 = s1 s1 = s pop = planet.pop t0 = planet.lastUpdate if s1 > pop : a0 = exp(rate*pop/s1) a1 = a0 + (t1 - t0) delta = s1*(log(a1/a0)/rate) if delta > s1 - pop : delta = s1 - pop planet.pop = pop + delta if pop > s1 : a0 = exp(rate*(1 - pop/s0)) a1 = a0 + (t1 - t0) delta = s0*(log(a1/a0)/rate) if delta > pop - s1 : delta = pop - s1 planet.pop = pop - delta planet.thresh = s1 planet.lastThresh = s0 planet.lastUpdate = t1

(les mises à jour sont bien sûr supposées en secondes)

Nimé, j'aimerais bien que tu le testes de ton côté pour des évolutions de seuils à ta convenance, et que tu compares un profil où un joueur fait beaucoup de mises à jour avec un où peu de mises à jour sont faites.

Nimeroni Mer 19 Juin - 17:45

Code:: public function calculSeuil(){ $seuil = $this->taille() ; if($this->consommationMedicament() == 1){ $joueur = $this->joueur() ; $bonusMedic = (float)(new Constante("planeteCroissanceMedicament"))->_valeur ; $medic = $this->ressource("medicament", $joueur) ; if($medic == NULL) $bonusMedic = 0 ; } else{ $bonusMedic = 0 ; } //!TODO: après la mise en place des races //if favorable(planet.climate, planet.owner) : //thresh = thresh * 2 $seuil = $seuil * (1 + $bonusMedic) ; return $seuil ; }

Code:: public function update($now = NULL) { if($now == NULL) $now = time() ; [...] //Calcul du temps $dernierUpdate = $this->dernierUpdate() ; $sec = $now - $dernierUpdate->format("U") ; [...] //Calcul de la population $taux = log((int) $tempsTaux->_valeur) ; //86400 $s0 = $this->dernierSeuil() ; $s1 = $this->seuilDemographique() ; $s = $this->calculSeuil() ; if($s != $s1){ $s0 = $s1 ; $s1 = $s ; } $pop = $this->population() ; if($s1 > $pop){ $a0 = exp($taux*$pop/$s1) ; $a1 = $a0 + $sec ; //$sec = t1 - t0 $delta = $s1*(log($a1/$a0)/$taux) ; if ($delta > ($s1 - $pop)) $delta = $s1 - $pop ; $this->population($pop + $delta) ; } if ($pop > $s1){ $a0 = exp($taux*(1 - $pop/$s0)) ; $a1 = $a0 + $sec ; $delta = $s0*(log($a1/$a0)/$taux) ; if ($delta > ($pop - $s1)) $delta = $pop - $s1 ; $this->population($pop + $delta) ; } $this->seuilDemographique(s1) ; $this->dernierSeuil(s0) ; $this->dernierUpdate(new DateTime(date("c", $now))) ; [...] }

Apeiron Jeu 20 Juin - 16:12

Pour la richesse le calcul est de la même forme, où le seuil est calculé par :

calculSeuilRichesse = 1 + (1 - taxation)/cteRich +50% si medocs, où :

- taxation est entre 0 et 1
- cteRich = 1 / (2*rmax/3 -1)
- la richesse maximale rmax est fixée pour le moment à 20.

Il faut rajouter à la structure de planète les variables suivantes (initialisées à la colo à) :

taxation (0)
richesse (1)
seuilRichesse (calculSeuilRichesse)
dernierSeuilRichesse (1)

Et le calcul est le même avec rich qui remplace pop, r0 qui remplace s0 et r1 qui remplace s1.
taux est le même, puisqu'on souhaite synchroniser les plafonds de la population et de la richesse.

Passons maintenant au calcul de l'imposition.
A un instant t, le gain est de taxation(t)*population(t)*richesse(t). Il faut donc calculer l'intégrale de cette quantité pour t parcourant la période de temps voulue.
taxation est considéré comme constant entre deux mises à jour, puisque changer la taxation met à jour la planète.

Pour la population, rappelons que :
(où pop0 est la population obtenue à la dernière mise à jour et time = 30)
- si s1 > pop0 alors pop(a) = s1*log(a)/taux si 1 < a < time, et s1 si a > time,
en partant de l'âge a0 = exp(taux*pop0/s1).
- si pop0 < s1 alors pop(a) = s0*(1 - log(a)/taux) si 1 < a < a2, et s1 si a > a2,
en partant de l'âge a0 = exp(taux*(1 - pop0/s0))
a2 est l'âge à partir duquel la population atteint s1, donné par :
a2 = exp(taux*(1 - s1/s0))

En partant d'une mise à jour t1 (en jours) et en sachant que la dernière mise à jour était t0, on a :

Code:: si s1 > pop0 a0 = exp(taux*pop0/s1) a1 = a0 + (t1 - t0) si a1 > time pop(a) = s1*log(a)/taux pour a allant de a0 à time pop(a) = s1 pour a allant de time à a1 sinon pop(a) = s1*log(a)/taux pour a allant de a0 à a1 si pop0 > s1 a0 = exp(taux*(1 - pop0/s0)) a1 = a0 + (t1 - t0) a2 = exp(taux*(1 - s1/s0)) si a1 > a2 pop(a) = s0*(1 - log(a)/taux) pour a allant de a0 à a2 pop(a) = s1 pour a allant de a2 à a1 sinon pop(a) = s0*(1 - log(a)/taux) pour a allant de a0 à a1

La richesse fonction selon les mêmes équations :
(où rich0 est la richesse obtenue à la dernière mise à jour, r0 l'ancien seuil de richesse et r1 le nouveau seuil de richesse)

Code:: si r1 > rich0 b0 = exp(taux*rich0/r1) b1 = b0 + (t1 - t0) si b1 > time rich(b) = r1*log(b)/taux pour b allant de b0 à time rich(b) = r1 pour b allant de time à b1 sinon rich(b) = r1*log(b)/taux pour b allant de b0 à b1 si rich0 > r1 b0 = exp(taux*(1 - rich0/r0)) b1 = b0 + (t1 - t0) b2 = exp(taux*(1 - r1/r0)) si b1 > b2 rich(b) = r0*(1 - log(b)/taux) pour b allant de b0 à b2 rich(b) = r1 pour b allant de b2 à b1 sinon rich(b) = r0*(1 - log(b)/taux) pour b allant de b0 à b1

Pour T = t1 - t0, on a log(a) = log(a0 + t) et log(b) = log(b0 + t) pour t entre 0 et T, donc :
imposition = integ[t = 0 à T]taxation(t)*population(t)*richesse(t)dt
              = taxation*integ[t = 0 à T]pop(a0 + t)*rich(b0 + t)dt

On se retrouve donc à intégrer un log(x)^2.
Après une courte révision de mes intégrations par parties, j'ai que :
- la primitive de log(x) est x*(log(x) -1)
- la primitive de (log(x))^2 est x*(log(x)^2 -2*log(x) +2)

Pour faire l'intégration, il faut de plus séparer les périodes logarithmiques des périodes constantes, en calculant a1 et b1 à time et a2 ou b2.
En tenant compte des différents cas possible, l'algo est donc :

imposition(planete, t1) =
   t0 = planete.derniereUpdate
   pop = planete.population          // population à t0
   s0 = planete.dernierSeuilPop
   s1 = planete.seuilPop
   rich = planete.richesse            // richesse à t0
   r0 = planete.dernierSeuilRich
   r1 = planete.seuilRich
   T = t1 - t0

   if s1 > pop
      a0 = exp(taux*pop/s1)
      if r1 > rich
         b0 = exp(taux*rich/r1)
         // la population et la richesse croissent
      if rich > r1
         b0 = exp(taux*(1 - rich/r0))
         b2 = exp(taux*(1 - r1/r0))
         // la population croît et la richesse décroît
   if pop > s1
      a0 = exp(taux*(1 - pop/s0))
      a2 = exp(taux*(1 - s1/s0))
      if r1 > rich
         b0 = exp(taux*rich/r1)
         // la population décroît et la richesse croît
      if rich > r1
         b0 = exp(taux*(1 - rich/r0))
         b2 = exp(taux*(1 - r1/r0))
         // la population et la richesse décroissent

Suite au prochain épisode...

Apeiron Dim 23 Juin - 18:06

Uniformisons un peu les formules...

Soit T = t1 - t0.
imposition = tax*intégrale[t = 0, T]pop(t)rich(t)dt
où pop(t) et rich(t) dépendent du cas :

si s1 > pop0
   a0 = exp(taux*pop0/s1)
   ta = time - a0
   pop(t) = s1*log(a0 + t)/taux pour t allant de 0 à min(T,ta)
si pop0 > s1
   a0 = exp(taux*(1 - pop0/s0))
   ta = exp(taux*(1 - s1/s0)) - a0
   pop(t) = s0*(1 - log(a0 + t)/taux) pour t allant de 0 à min(T,ta)
dans les deux cas, si T > ta, il faut rajouter :
   pop(t) = s1 pour t allant de ta à T
enfin, si pop0 = s1
   pop(t) = s1 pour t allant de 0 à T

si r1 > rich0
   b0 = exp(taux*rich0/r1)
   tb = time - b0
   rich(t) = r1*log(b0 + t)/taux pour t allant de 0 à min(T,tb)
si rich0 > r1
   b0 = exp(taux*(1 - rich0/r0))
   tb = exp(taux*(1 - r1/r0)) - b0
   rich(t) = r0*(1 - log(b0 + t)/taux) pour t allant de 0 à min(T,tb)
dans les deux cas, si T > tb, il faut rajouter :
   rich(t) = r1 pour t allant de tb à T
enfin, si rich0 = r1
   rich(t) = r1 pour t allant de 0 à T

J'obtiens donc l'algo suivant :

Code:: // cette initialisation sert en cas de non évolution entre t0 et t1 ta = 0 tb = 0 // les morceaux sont ajoutés au fur et à mesure, tax étant appliqué à la fin imposition = 0 // préparation des constantes selon les cas si s1 > pop0 a0 = exp(taux*pop0/s1) ta = time - a0 si pop0 > s1 a0 = exp(taux*(1 - pop0/s0)) ta = exp(taux*(1 - s1/s0)) - a0 si r1 > rich0 b0 = exp(taux*rich0/r1) tb = time - b0 si rich0 > r1 b0 = exp(taux*(1 - rich0/r0)) tb = exp(taux*(1 - r1/r0)) - b0 maxt = max(ta,tb) mint = min(ta,tb) // cas sans log si maxt < T imposition = imposition + s1*r1*(T - maxt) // cas avec un seul log si mint < maxt si mint = ta integ = (b0 + maxt)*log(b0 + maxt) - (b0 + mint)*log(b0 + mint) - (maxt - mint) si r1 > rich0 imposition = imposition + s1*r1*integ/taux si rich0 > r1 imposition = imposition + s1*r0*((maxt - mint) - integ/taux) si mint = tb integ = (a0 + maxt)*log(a0 + maxt) - (a0 + mint)*log(a0 + mint) - (maxt - mint) si s1 > pop0 imposition = imposition + r1*s1*integ/taux si pop0 > s1 imposition = imposition + r1*s0*((maxt - mint) - integ/taux) // cas avec deux logs si 0 < mint intega = ((a0 + mint)*log(a0 + mint) - a0*log(a0) - mint)/taux integb = ((b0 + mint)*log(b0 + mint) - b0*log(b0) - mint)/taux integab = (((a0 + b0)/2 + mint)*log(a0 + mint)*log(b0 + mint) - (a0 + mint)*log(a0 + mint) - (b0 + mint)*log(b0 + mint)+ 2*mint - ((a0 + b0)/2)*log(a0)*log(b0) + a0*log(a0) + b0*log(b0))/(taux^2) si s1 > pop0 si r1 > rich0 imposition = imposition + s1*r1*integab si rich0 > r1 imposition = imposition + s1*r0*(intega - integab) si pop0 > s1 si r1 > rich0 imposition = imposition + s0*r1*(integb - integab) si rich0 > r1 imposition = imposition + s0*r0*(1 - intega - integb + integab) imposition = tax*imposition

Le post suivant explique l'origine de cette expression pour integab. Les tests numériques seront faits plus tard.

Apeiron Lun 24 Juin - 13:43

integab = integ[t = 0 à mint]log(a0 + t)*log(b0 + t)/(taux^2) dt

La primitive de log[a + t]*log[b + t] est :
((a + t)*log[a + t] - t)*(log[b + t] - 1) + t - b*log[b + t] + (b - a)*prim[log[a + t]/(b + t)]

Quelle est la primitive de log[a + t]/(b + t) ?
D'après Wolgram Mathematica la réponse est : Log[a + t]*Log[1 - (a + t)/(a - b)] + diLog[(a + t)/(a - b)].

Moi, quand je prends comme axiome que la primitive de log[x]/(x - 1) est - diLog[1 - x] je tombe sur une primitive de : log[a - b]*log[(a + t)/(a - b) -1] - diLog[1 - (a + t)/(a - b)], ce qui m'embête.

Je ne crois pas sur parole Wolgram parce que je l'ai déjà vu faire apparaître un a que je n'explique pas, mais là c'est un peu gros, à moins qu'il y ait des relations avec le diLog que j'ignore.

Pour l'approximation du diLog dans les expressions diLog[(a + t)/(a - b)] et diLog[1 - (a + t)/(a - b)] :
Encadrons tout ça : 1 < a,b < 30 et 0 < t < 2 (car t est en jours, et on impose une mise à jour au moins toutes les 48h).
On peut supposer a != b car sinon l'intégrale est facile à traiter [1]. On a donc |a - b| > p, où p est la précision du test d'égalité.
Donc :
- si a > b alors 1/29 < (a + t)/(a - b) < 32/p et (p - 32)/p < 1 - (a + t)/(a - b) < 28/29
- si a < b alors -1/p < (a + t)/(a - b) < -32/29 et 61/29 < 1 - (a + t)/(a - b) < (1+p)/p
L'encadrement n'est pas génial...

Quand on voit la tête de prim[log[a + t]/(b + t)], avec des logs et diLog en sens inverses, il est permis de se demander s'il n'y aurait pas des simplifications non pas formelles mais numériques. Dit autrement : est-ce que cette quantité est approximable par une constante (par exemple 0) ? Si c'était le cas, on n'aurait qu'à la supprimer de l'équation et dire aux joueurs qu'il y a un peu d'aléatoire (en fait une absence de bruit Very Happy

) dans l'imposition. Mais bon, pour aller plus loin dans cette voie il me faudrait un logiciel approprié...

[1] Car dans ce cas on aurait :
integab = integ[t = 0 à maxt]log(a + t)*log(b + t)/(taux^2) dt = integ[t = 0 à maxt]log(a + t)^2/(taux^2) dt
Or, la primitive de log(a + t)^2 est : (a + t)*((Log[a + t] - 2)*Log[a + t] + 2), soit en développant :
(a + t)*log[a + t]^2 - 2*(a + t)*log[a + t] + 2*(a + t)
ce qui est très proche de la primitive de log[a + t]*log[b + t] où on a fait a = b :
(a + t)*log[a + t]^2 - 2*(a + t)*log[a + t] + 2*t
en négligeant (b - a)*prim[log[a + t]/(b + t)].
(et encore plus si on croit Wolgram parce qu'il fait apparaître un +a)

Du coup, même sans logiciel, j'ai bien envie d'approximer la primitive de log[a + t]*log[b + t] par :
(a + t)*log[a + t]*log[b + t] - (a + t)*log[a + t] - (b + t)*log[b + t]+ 2*t + a
Ou peut-être (pour la symétrie en a et b) :
((a + b)/2 + t)*log[a + t]*log[b + t] - (a + t)*log[a + t] - (b + t)*log[b + t]+ 2*t + a + b

PS :
C'est une primitive donc les constantes sont à supprimer de toute façon. Ensuite, pour le log[a + t]*log[b + t], on cherche à voir comment moyenniser le a avec le b.
(a + t)*log[a + t]*log[b + t] + (b + t)*log[a + t]*log[b + t] = (a + b + 2*t)*log[a + t]*log[b + t] donc prendre (a + b)/2 est correct (on aurait aussi pu prendre la racine de ab, mais ici la moyenne arithmétique colle mieux que la moyenne géométrique).

Conclusion :
Je vais prendre comme primitive de log(a0 + t)*log(b0 + t) :
((a0 + b0)/2 + t)*log(a0 + t)*log(b0 + t) - (a0 + t)*log(a0 + t) - (b0 + t)*log(b0 + t)+ 2*t
donc :
integab = ((a0 + b0)/2 + mint)*log(a0 + mint)*log(b0 + mint) - (a0 + mint)*log(a0 + mint) - (b0 + mint)*log(b0 + mint)+ 2*mint - ((a0 + b0)/2)*log(a0)*log(b0) + a0*log(a0) + b0*log(b0)

(le post précédent a donc été édité avec cette expression, les essais numériques suivront plus tard)

Apeiron Sam 29 Juin - 21:14

Voici le code que je teste actuellement :

Code:: //définition des constantes et des éléments pertinents d'une planète float time = 30.0; float taux = log(30.0); typedef struct planete planete; struct planete { float population; float seuilPop; float dernierSeuilPop; float richesse; float seuilRich; float dernierSeuilRich; float taxation; float derniereUpdate; }; //Nous pourrons faire varier le seuil et la richesse en fonction de la date float calculSeuilPop(planete p, float now){ float seuilPop = 9.0; if(now > 30.0){seuilPop = 0.0;} return seuilPop; } float calculSeuilRich(planete p, float now){ float seuilRich = 21.0; if(now > 30.0){seuilRich = 1.0;} return seuilRich; } //colo initialise les éléments d'une planète // posons pour l'instant uen taxation de 50% planete colo(planete p, float now){ p.population = 0.0; p.seuilPop = calculSeuilPop(p,now); p.dernierSeuilPop = 0.0; p.richesse = 1.0; p.seuilRich = calculSeuilRich(p,now); p.dernierSeuilRich = 1.0; p.taxation = 0.5; p.derniereUpdate = now; return p; } planete updatePop(planete p, float now){ float pop = p.population; float s0 = p.dernierSeuilPop; float s1 = p.seuilPop; float t0 = p.derniereUpdate; float t1 = now; //mise à jour éventuelle du seuil float s = calculSeuilPop(p,now); if(s != s1){ s0 = s1; s1 = s ; } float a0, a1, delta; // croissance de la population if(s1 > pop){ a0 = exp(taux*pop/s1); a1 = a0 + (t1 - t0); delta = s1*(log(a1/a0)/taux); if(delta > (s1 - pop)){delta = s1 - pop;} pop = pop + delta; } // décroissance de la population if(pop > s1){ a0 = exp(taux*(1 - pop/s0)); a1 = a0 + (t1 - t0); delta = s0*(log(a1/a0)/taux); if(delta > (pop - s1)){delta = pop - s1;} pop = pop - delta; } // mise à jour de l'état p.population = pop; p.seuilPop = s1; p.dernierSeuilPop = s0; return p; } planete updateRich(planete p, float now){ float rich = p.richesse; float r0 = p.dernierSeuilRich; float r1 = p.seuilRich; float t0 = p.derniereUpdate; float t1 = now; //mise à jour éventuelle du seuil float r = calculSeuilRich(p,now); if(r != r1){ r0 = r1; r1 = r ; } float b0, b1, delta; // croissance de la richesse if(r1 > rich){ b0 = exp(taux*rich/r1); b1 = b0 + (t1 - t0); delta = r1*(log(b1/b0)/taux); if(delta > (r1 - rich)){delta = r1 - rich;} rich = rich + delta; } // décroissance de la richesse if(rich > r1){ b0 = exp(taux*(1 - rich/r0)); b1 = b0 + (t1 - t0); delta = r0*(log(b1/b0)/taux); if(delta > (rich - r1)){delta = rich - r1;} rich = rich - delta; } // mise à jour de l'état p.richesse = rich; p.seuilRich = r1; p.dernierSeuilRich = r0; return p; } float imposition(planete p, float now){ float t0 = p.derniereUpdate; float t1 = now; float T = t1 - t0; // cette initialisation sert en cas de non évolution entre t0 et t1 float ta = 0.0; float tb = 0.0; // les morceaux sont ajoutés au fur et à mesure, la taxation étant appliquée à la fin float impot = 0.0; float pop0 = p.population; float rich0 = p.richesse; //cette partie-là mériterait une optimisation, vu qu'elle copie les mises à jours de la pop et de la richesse float s0 = p.dernierSeuilPop; float s1 = p.seuilPop; float s = calculSeuilPop(p,now); if(s != s1){ s0 = s1; s1 = s ; } float r0 = p.dernierSeuilRich; float r1 = p.seuilRich; float r = calculSeuilRich(p,now); if(r != r1){ r0 = r1; r1 = r ; } // préparation des constantes selon les cas float a0, b0; if(s1 > pop0){ a0 = exp(taux*pop0/s1); ta = time - a0; } if(pop0 > s1){ a0 = exp(taux*(1.0 - pop0/s0)); ta = exp(taux*(1.0 - s1/s0)) - a0; } if(r1 > rich0){ b0 = exp(taux*rich0/r1); tb = time - b0; } if(rich0 > r1){ b0 = exp(taux*(1.0 - rich0/r0)); tb = exp(taux*(1.0 - r1/r0)) - b0; } float maxt = fmax(ta,tb); float mint = fmin(ta,tb); // cas sans log if(maxt < T){ impot = impot + s1*r1*(T - maxt); } // cas avec un seul log float integ; if(mint < maxt){ if(mint = ta){ integ = (b0 + maxt)*log(b0 + maxt) - (b0 + mint)*log(b0 + mint) - (maxt - mint); if(r1 > rich0){ impot = impot + s1*r1*integ/taux; } if(rich0 > r1){ impot = impot + s1*r0*((maxt - mint) - integ/taux); } } if(mint = tb){ integ = (a0 + maxt)*log(a0 + maxt) - (a0 + mint)*log(a0 + mint) - (maxt - mint); if(s1 > pop0){ impot = impot + r1*s1*integ/taux; } if(pop0 > s1){ impot = impot + r1*s0*((maxt - mint) - integ/taux); } } } // cas avec deux logs float intega, integb, integab; if(0.0 < mint){ intega = ((a0 + mint)*log(a0 + mint) - a0*log(a0) - mint)/taux; integb = ((b0 + mint)*log(b0 + mint) - b0*log(b0) - mint)/taux; integab = (((a0 + b0)/2.0 + mint)*log(a0 + mint)*log(b0 + mint) - (a0 + mint)*log(a0 + mint) - (b0 + mint)*log(b0 + mint)+ 2.0*mint - ((a0 + b0)/2.0)*log(a0)*log(b0) + a0*log(a0) + b0*log(b0))/pow(taux,2); if(s1 > pop0){ if(r1 > rich0){ impot = impot + s1*r1*integab; } if(rich0 > r1){ impot = impot + s1*r0*(intega - integab); } } if(pop0 > s1){ if(r1 > rich0){ impot = impot + s0*r1*(integb - integab); } if(rich0 > r1){ impot = impot + s0*r0*(1.0 - intega - integb + integab); } } } impot = p.taxation*impot; return impot; } int main(){ //le test se déroule entre tmin et tmax float tmin = 0.0; float tmax = 60.0; int nb_updates = 10; //initialisation de la planète à la date 0 float pop, rich; float impot = 0.0; planete p = colo(p,0.0); //t est la date de chaque mise à jour float t; float pas = (tmax - tmin)/nb_updates; for(t = tmin + pas ; t <= tmax ; t = t + pas){ impot = impot + imposition(p,t); p = updatePop (p,t); p = updateRich(p,t); //mémoriser la derniere update doit se faire ici p.derniereUpdate = t; /* test au fil du temps de la population et de la richesse pop = p.population; rich = p.richesse; printf("date : %f\n",t); printf("pop : %f\n",pop); printf("rich : %f\n\n",rich); */ } printf("updates : %d\n",nb_updates); printf("impot : %f\n\n",impot); return 0; }

L'évolution de la population ainsi que de la richesse est bonne, mais le calcul de l'imposition est complètement aberrant.

Apeiron Lun 8 Juil - 1:32

Pendant que je serai à Kandorya, tu peux implémenter l'évolution de la population et de la richesse, ainsi que l'interface de taxation et l'affichage prédictif. L'imposition devra attendre que je revienne et que je résolve le problème. Après on passera à la production et la consommation de ressources, voire au crime.

Contenu sponsorisé

Quelques calculs...

Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...

Re: Quelques calculs...