Add absolute tolerance to KMeans algorithm

IBPSA/Fluid/Geothermal/Borefields/BaseClasses/HeatTransfer/ThermalResponseFactors/Validation/ClusterBoreholes_100boreholes.mo

@@ -15,10 +15,11 @@ model ClusterBoreholes_100boreholes
   parameter Integer cluSiz[k](each fixed=false) "Size of the clusters";
 
   parameter Integer labelsExp[nBor]=
-    {3, 3, 3, 4, 4, 4, 4, 3, 3, 3, 3, 4, 4, 2, 2, 2, 2, 4, 4, 3, 3, 4, 2, 2, 1,
-     1, 2, 2, 4, 3, 4, 2, 2, 1, 1, 1, 1, 2, 2, 4, 4, 2, 1, 1, 1, 1, 1, 1, 2, 4,
-     4, 2, 1, 1, 1, 1, 1, 1, 2, 4, 4, 2, 2, 1, 1, 1, 1, 2, 2, 4, 3, 4, 2, 2, 1,
-     1, 2, 2, 4, 3, 3, 4, 4, 2, 2, 2, 2, 4, 4, 3, 3, 3, 3, 4, 4, 4, 4, 3, 3, 3}
+    {1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 2, 4, 4, 4, 4, 4, 4, 2, 1, 2, 4, 4, 3,
+     3, 3, 3, 4, 4, 2, 2, 4, 3, 3, 3, 3, 3, 3, 4, 2, 2, 4, 3, 3, 3, 3, 3, 3,
+     4, 2, 2, 4, 3, 3, 3, 3, 3, 3, 4, 2, 2, 4, 3, 3, 3, 3, 3, 3, 4, 2, 2, 4,
+     4, 3, 3, 3, 3, 4, 4, 2, 1, 2, 4, 4, 4, 4, 4, 4, 2, 1, 1, 1, 2, 2, 2, 2,
+     2, 2, 1, 1}
     "Expected cluster labels";
 
   // Comparison result

IBPSA/Resources/ReferenceResults/Dymola/IBPSA_Fluid_Geothermal_Borefields_BaseClasses_HeatTransfer_ThermalResponseFactors_Validation_GFunction_100boreholes.txt

@@ -1,9 +1,9 @@
-last-generated=2022-09-07
+last-generated=2025-03-29
 statistics-simulation=
 {
   "linear": " ",
   "nonlinear": " ",
   "numerical Jacobians": "0"
 }
 time=[0e+00, 3.70000003072e+11]
-gFun_int=[0e+00, 5.314130027278435e+01, 5.783808902952384e+01, 5.949810404381666e+01, 6.029800416918814e+01, 6.07697754272593e+01, 6.105554196943749e+01, 6.12624855459325e+01, 6.13997192460549e+01, 6.151123426988645e+01, 6.159519197084513e+01, 6.165727618463693e+01, 6.17088775545948e+01, 6.175517273659016e+01, 6.179167176974389e+01, 6.182075881365947e+01, 6.184498977959359e+01, 6.186692811114929e+01, 6.188775634245692e+01, 6.190541076800866e+01, 6.192027664779165e+01, 6.193291092840229e+01, 6.194388200029567e+01, 6.195374678888403e+01, 6.196306609780395e+01, 6.197225189208984e+01, 6.198055649077359e+01, 6.19878959704005e+01, 6.199438858672767e+01, 6.200014876881503e+01, 6.200530242701172e+01, 6.200997161825369e+01, 6.201427459826257e+01, 6.201832962282079e+01, 6.202226257705131e+01, 6.202618790189738e+01, 6.202994918630674e+01, 6.203339004457818e+01, 6.203654480033171e+01, 6.203942871239696e+01, 6.204207611306958e+01, 6.204450989058045e+01, 6.204675674782058e+01, 6.204883957299383e+01, 6.205078888369859e+01, 6.205262374560362e+01, 6.205437469236476e+01, 6.205606460392911e+01, 6.205772018316401e+01, 6.205936813296384e+01, 6.20610237121582e+01, 6.206260681205722e+01, 6.206409072977328e+01, 6.20654869093779e+01, 6.20667953509319e+01, 6.206802749843249e+01, 6.206918716665799e+01, 6.20702743556287e+01, 6.207130431967757e+01, 6.20722732532121e+01, 6.207318878094436e+01, 6.207405853225876e+01, 6.207488632187959e+01, 6.207567977919125e+01, 6.207644271890113e+01, 6.207717895572002e+01, 6.207789611906271e+01, 6.20786018382528e+01, 6.207929611339839e+01, 6.20799865738467e+01, 6.208068084899229e+01, 6.208137512109734e+01, 6.208203887870534e+01, 6.208266830399742e+01, 6.208326339700734e+01, 6.208382415771484e+01, 6.208435821551446e+01, 6.208486175572244e+01, 6.208534240771305e+01, 6.208579635682279e+01, 6.20862312324327e+01, 6.208664321984549e+01, 6.208703613375845e+01, 6.208741378892288e+01, 6.20877799998739e+01, 6.208813095214395e+01, 6.208847046032221e+01, 6.208880233902487e+01, 6.208913040310457e+01, 6.208945083778975e+01, 6.20897674549736e+01, 6.209008788959122e+01, 6.209040450963318e+01, 6.209072875896157e+01, 6.209105682300074e+01, 6.209138870181825e+01, 6.209173202466e+01, 6.209208679162732e+01, 6.209245681735668e+01, 6.209283828721838e+01, 6.209323883056641e+01]
+gFun_int=[0e+00, 5.325750358097349e+01, 5.795590595455504e+01, 5.961476130093882e+01, 6.041376878708399e+01, 6.088492587886667e+01, 6.117026136034042e+01, 6.137689213155998e+01, 6.151390076459565e+01, 6.162522886829792e+01, 6.170904542540963e+01, 6.177102282763745e+01, 6.182252883024646e+01, 6.186874390356883e+01, 6.190517808684201e+01, 6.193421172504487e+01, 6.195840072929303e+01, 6.198029709915853e+01, 6.200108718353151e+01, 6.201870727679434e+01, 6.203354263899414e+01, 6.20461578460965e+01, 6.205710602981472e+01, 6.206695174495223e+01, 6.207625579506957e+01, 6.208541870117188e+01, 6.209371185575707e+01, 6.210103607660167e+01, 6.210751343413978e+01, 6.211326598683261e+01, 6.211841201563477e+01, 6.212306976278494e+01, 6.212736129870709e+01, 6.213141250855115e+01, 6.213533783338038e+01, 6.213925552885051e+01, 6.214300918386534e+01, 6.214644241274562e+01, 6.214958953910124e+01, 6.215246963646923e+01, 6.215511322244458e+01, 6.215754318523454e+01, 6.215978241310378e+01, 6.216186142357977e+01, 6.216380691958726e+01, 6.216563796679502e+01, 6.216738509888592e+01, 6.216907119574625e+01, 6.217072677494737e+01, 6.217237091006345e+01, 6.217402267456055e+01, 6.21756019597623e+01, 6.217708206278109e+01, 6.217847442770872e+01, 6.217978286921542e+01, 6.218101501674303e+01, 6.218217087027127e+01, 6.218325805924198e+01, 6.218428039392673e+01, 6.218524932743085e+01, 6.218616485516311e+01, 6.218703079179038e+01, 6.218785858140107e+01, 6.218865203870936e+01, 6.218941116372535e+01, 6.219014740054423e+01, 6.219086456386328e+01, 6.219156646837975e+01, 6.219225692886523e+01, 6.219294738927639e+01, 6.219363784977539e+01, 6.219433212182976e+01, 6.219499587943776e+01, 6.219562149004609e+01, 6.219621658304925e+01, 6.219677734375e+01, 6.219731140154962e+01, 6.21978149417576e+01, 6.219829177907121e+01, 6.219874572816068e+01, 6.219918060377059e+01, 6.219959259118338e+01, 6.219998550509634e+01, 6.220036316020671e+01, 6.220072555651452e+01, 6.220107650878457e+01, 6.220141601688852e+01, 6.220174789566549e+01, 6.220207214504793e+01, 6.220239257973311e+01, 6.220271301154666e+01, 6.220302963153458e+01, 6.220334625152249e+01, 6.220366668625496e+01, 6.220399475028737e+01, 6.220433044372783e+01, 6.220467376657633e+01, 6.220502471885315e+01, 6.220539092991903e+01, 6.220577239976721e+01, 6.220617294311523e+01]

IBPSA/Resources/ReferenceResults/Dymola/IBPSA_Fluid_Geothermal_Borefields_Validation_ConstantHeatInjection_100Boreholes.txt

@@ -1,10 +1,10 @@
-last-generated=2019-05-31
+last-generated=2025-03-29
 statistics-simulation=
 {
   "linear": "0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0",
   "nonlinear": " ",
   "number of continuous time states": "42",
   "numerical Jacobians": "0"
 }
-borHol.TBorAve=[2.731499938964844e+02, 2.969233097411025e+02, 3.061154789798177e+02, 3.117919010735677e+02, 3.156233524842303e+02, 3.185083012144097e+02, 3.206774600774739e+02, 3.223859561644965e+02, 3.237309573611698e+02, 3.249128115592614e+02, 3.258684390314781e+02, 3.266721499351128e+02, 3.273854677675781e+02, 3.279972231598931e+02, 3.285019228322349e+02, 3.289506535865341e+02, 3.293434755466944e+02, 3.297362669857326e+02, 3.300673524590135e+02, 3.3034216348059e+02, 3.306169435382008e+02, 3.308567504599264e+02, 3.310732423402778e+02, 3.312897037022963e+02, 3.314708253237771e+02, 3.316351623535156e+02, 3.317994996471114e+02, 3.31935394501481e+02, 3.320621033808788e+02, 3.321888122720659e+02, 3.322926638380795e+02, 3.323919984610638e+02, 3.324913025909467e+02, 3.325905761996577e+02, 3.326898803225085e+02, 3.327891847630003e+02, 3.328550720996233e+02, 3.329205932864627e+02, 3.329860841751099e+02, 3.33051574845879e+02, 3.331170655130119e+02, 3.331804199496577e+02, 3.332311401228471e+02, 3.332818299520609e+02, 3.333325196012856e+02, 3.333832092615076e+02, 3.334338990923591e+02, 3.334806214316239e+02, 3.335180664687611e+02, 3.335555115023921e+02, 3.335929260253906e+02, 3.336303711927259e+02, 3.336678162241643e+02, 3.337008667235286e+02, 3.337286683031417e+02, 3.337564391916251e+02, 3.33784240771276e+02, 3.338120423508504e+02, 3.338398132393775e+02, 3.338638611447303e+02, 3.338855286755524e+02, 3.339071655585993e+02, 3.339288330806903e+02, 3.339504699628948e+02, 3.339721374974545e+02, 3.339937745112744e+02, 3.340154419049498e+02, 3.340370789273128e+02, 3.34058746313598e+02, 3.340803833326005e+02, 3.341020508645744e+02, 3.341164245588113e+02, 3.341301575019587e+02, 3.341438904450156e+02, 3.341576538190172e+02, 3.341713867621481e+02, 3.341851197052002e+02, 3.341988525616359e+02, 3.342125855047053e+02, 3.342263488786904e+02, 3.342400818218279e+02, 3.342538147648898e+02, 3.342668151959654e+02, 3.342773132792386e+02, 3.342878112723611e+02, 3.34298309347003e+02, 3.343088074320851e+02, 3.343193054251311e+02, 3.34329833984375e+02, 3.343403320786896e+02, 3.343508300920124e+02, 3.343613281501781e+02, 3.343718262343417e+02, 3.343823242483112e+02, 3.343912964283731e+02, 3.343989563040388e+02, 3.344066162317745e+02, 3.344142761603534e+02, 3.344219360351562e+02, 3.344295959831346e+02, 3.34437255859375e+02]
 time=[0e+00, 3.1536001024e+10]
+borHol.TBorAve=[2.731499938964844e+02, 2.970736698528646e+02, 3.063688664364149e+02, 3.121055912630208e+02, 3.159714359837963e+02, 3.188781437463107e+02, 3.210608218951823e+02, 3.227778629050565e+02, 3.241284488197207e+02, 3.253146975468917e+02, 3.262730105646812e+02, 3.27078674593316e+02, 3.277935183046875e+02, 3.284063418122369e+02, 3.289115908009849e+02, 3.29360779318956e+02, 3.297539674900538e+02, 3.301471251434999e+02, 3.304784547539354e+02, 3.307533573218809e+02, 3.310282289385914e+02, 3.312680358608957e+02, 3.314845277430555e+02, 3.31701019620265e+02, 3.318820802065896e+02, 3.320463256835938e+02, 3.322105714244551e+02, 3.323463442085122e+02, 3.324729310175975e+02, 3.325994873208941e+02, 3.327032778448076e+02, 3.328024904165695e+02, 3.329017029850446e+02, 3.330009155516536e+02, 3.331001281217726e+02, 3.331993410130003e+02, 3.332651673066532e+02, 3.333305969485721e+02, 3.333959962844849e+02, 3.334614259200977e+02, 3.335268555555463e+02, 3.33590148953564e+02, 3.336408080915971e+02, 3.336914368856547e+02, 3.337420654997231e+02, 3.337926941247889e+02, 3.3384335341724e+02, 3.338900147343597e+02, 3.339273987441517e+02, 3.339647827426265e+02, 3.340021667480469e+02, 3.340395508802259e+02, 3.340769348765081e+02, 3.341099548582943e+02, 3.34137695402751e+02, 3.341654662894986e+02, 3.341932068357291e+02, 3.342209473801473e+02, 3.34248687754574e+02, 3.342727051388709e+02, 3.342943421521149e+02, 3.343159485097698e+02, 3.343375855220966e+02, 3.343591918867229e+02, 3.343808288944479e+02, 3.344024353999462e+02, 3.344240417619382e+02, 3.34445678770844e+02, 3.344672851354277e+02, 3.344889221509599e+02, 3.345105286477776e+02, 3.345248718244363e+02, 3.345386047675837e+02, 3.345523071930625e+02, 3.345660400494859e+02, 3.345797424750388e+02, 3.345934754180908e+02, 3.346071777569484e+02, 3.346209107000178e+02, 3.346346130388466e+02, 3.346483459819842e+02, 3.346620484074679e+02, 3.346750183209654e+02, 3.346855163574219e+02, 3.346959838832508e+02, 3.347064819544249e+02, 3.347169495219289e+02, 3.347274475136726e+02, 3.347379150772095e+02, 3.347484131333771e+02, 3.347588806291218e+02, 3.347693787124655e+02, 3.347798767714511e+02, 3.347903442678424e+02, 3.347992859303262e+02, 3.34806945805992e+02, 3.348145752369866e+02, 3.348222351447284e+02, 3.348298950195312e+02, 3.348375244499315e+02, 3.348451843261719e+02]

IBPSA/Utilities/Clustering/KMeans.mo

@@ -6,6 +6,7 @@ impure function KMeans "k-means clustering algorithm"
   input Integer n_samples "Number of samples";
   input Integer n_features "Number of features";
   input Real relTol=1e-5 "Relative tolerance on cluster positions";
+  input Real absTol=1e-8 "Absolute tolerance on cluster positions";
   input Integer max_iter=500 "Maximum number of k-means iterations";
   input Integer n_init=10 "Number of runs with randomized centroid seeds";
   input Integer n_cluster_size=0 "Length of the cluster_size output vector";
@@ -16,7 +17,8 @@ impure function KMeans "k-means clustering algorithm"
 protected
   Real old_centroids[n_clusters,n_features] "Previous iteration centroids";
   Real new_centroids[n_clusters,n_features] "Next iteration centroids";
-  Real delta_centroids "Maximum relative displacement of cluster centroids between two k-means iterations";
+  Real relDelta_centroids "Maximum relative displacement of cluster centroids between two k-means iterations";
+  Real absDelta_centroids "Maximum absolute displacement of cluster centroids between two k-means iterations";
   Integer new_labels[n_samples] "Next iteration cluster labels";
   Real new_inertia "Inertia of the samples during the current run";
   Real inertia "Minimum inertia of the samples since first run";
@@ -31,30 +33,30 @@ algorithm
  id := Modelica.Math.Random.Utilities.initializeImpureRandom(seed);
 
   // ---- Perform n_init successive runs of the k-means algorithm
- for run in 1:n_init loop
+  inertia := Modelica.Constants.inf;
+  for run in 1:n_init loop
     // ---- Select initial centroids at random
     // Select 3 non-repeated data points in the data set
     n := Modelica.Math.Random.Utilities.impureRandomInteger(id,1,n_samples);
     old_centroids[1,:] := data[n,:];
     for i in 2:n_clusters loop
-      n := Modelica.Math.Random.Utilities.impureRandomInteger(id,1,n_samples);
-      old_centroids[i,:] := data[n,:];
-      min_dis := Modelica.Math.Vectors.norm(old_centroids[i,:]-old_centroids[1,:], p=2)^2;
+      min_dis := 0;
       while min_dis < Modelica.Constants.eps loop
         n := Modelica.Math.Random.Utilities.impureRandomInteger(id,1,n_samples);
         old_centroids[i,:] := data[n,:];
         min_dis := Modelica.Math.Vectors.norm(old_centroids[i,:]-old_centroids[1,:], p=2)^2;
-        for j in 1:i-1 loop
-          dis := Modelica.Math.Vectors.norm(old_centroids[j,:]-old_centroids[i,:], p=2)^2;
+        for j in 2:i-1 loop
+          dis := Modelica.Math.Vectors.norm(old_centroids[i,:]-old_centroids[j,:], p=2)^2;
           min_dis := min(dis, min_dis);
         end for;
       end while;
     end for;
 
     // ---- k-means iterations
     k_iter := 0;
-    delta_centroids := 2*relTol;
-    while k_iter < max_iter and delta_centroids > relTol loop
+    relDelta_centroids := 2*relTol;
+    absDelta_centroids := 2*absTol;
+    while k_iter < max_iter and (relDelta_centroids > relTol or absDelta_centroids > absTol) loop
       k_iter := k_iter + 1;
 
       // Find centroid closest to each data point
@@ -71,7 +73,8 @@ algorithm
       end for;
 
       // Re-evaluate position of the centroids
-      delta_centroids := 0;
+      relDelta_centroids := 0;
+      absDelta_centroids := 0;
       for j in 1:n_clusters loop
         n := sum(if new_labels[i]==j then 1 else 0 for i in 1:n_samples);
         new_centroids[j,:] := zeros(n_features);
@@ -84,7 +87,8 @@ algorithm
         else
           new_centroids[j,:] := old_centroids[j,:];
         end if;
-        delta_centroids := max(delta_centroids, sum((new_centroids[j,:] - old_centroids[j,:])./old_centroids[j,:]));
+        relDelta_centroids := max(relDelta_centroids, sum(abs(new_centroids[j,:] - old_centroids[j,:]) ./ (abs(old_centroids[j,:]) .+ Modelica.Constants.eps)));
+        absDelta_centroids := max(absDelta_centroids, sum(abs(new_centroids[j,:] - old_centroids[j,:])));
       end for;
       old_centroids := new_centroids;
     end while;
@@ -93,7 +97,7 @@ algorithm
     new_inertia := 0;
     for i in 1:n_samples loop
       dis := Modelica.Math.Vectors.norm(data[i,:]-centroids[new_labels[i],:], p=2)^2;
-      new_inertia := inertia + dis;
+      new_inertia := new_inertia + dis;
     end for;
 
     // Keep run results if inertia is minimum
@@ -129,6 +133,15 @@ modifying the constant <code>seed</code>.
 </html>", revisions="<html>
 <ul>
 <li>
+March 18, 2025 by Massimo Cimmino<br/>
+Added absolute tolerance. The algorithm stops when any of the relative and
+absolute tolerances is satisfied. This fixes errors that occur when a centroid
+has a value close to zero on any of its axes. See
+<a href=\"https://github.com/ibpsa/modelica-ibpsa/issues/1985\">#1985</a>.
+Fixed the initial selection of centroids to avoid repeated centroids. See also
+<a href=\"https://github.com/ibpsa/modelica-ibpsa/issues/1976\">#1976</a>.
+</li>
+<li>
 February 1, 2023, by Michael Wetter:<br/>
 Added <code>impure</code> declaration which is needed for compliance with the Modelica Language Specification,
 and is required by Optimica.

IBPSA/Utilities/Clustering/Validation/KMeans_1d.mo

@@ -15,7 +15,7 @@ model KMeans_1d "Model that verifies the k-means clustering function for 1d data
     "Size of the clusters";
 
   parameter Integer labelsExp[nDat]=
-    {2,2,2,1,3}
+    {3,2,2,1,1}
     "Expected cluster labels";
 
   // Comparison result

IBPSA/Utilities/Clustering/Validation/KMeans_2d.mo

@@ -22,7 +22,7 @@ model KMeans_2d "Model that verifies the k-means clustering function for 2d data
     "Size of the clusters";
 
   parameter Integer labelsExp[nDat]=
-    {2,2,1,2,1,3}
+    {3,3,2,2,1,1}
     "Expected cluster labels";
 
   // Comparison result